1 DATA

1.1 Pop sive over time dataset

data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
  
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5
## 1          no
## 2          no
## 3          no
## 4          no
## 5          no
## 6          no
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571               Med                  5               6     8/25/21
## 572               Med                  5               6     8/25/21
## 573               Med                  5               6     8/25/21
## 574               Med                  5               6     8/25/21
## 575               Med                  5               6     8/25/21
## 576               Med                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5
## 571          no                     no         yes
## 572     extinct                    yes         yes
## 573          no                     no          no
## 574     extinct                    yes         yes
## 575          no                     no          no
## 576          no                     no          no
dim(data)
## [1] 576  23
summary(data)
##     Block            Old_Label            Label            Population       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  Genetic_Diversity  Generation_Parents Generation_Eggs Date_Census       
##  Length:576         Min.   :0.0        Min.   :1.0     Length:576        
##  Class :character   1st Qu.:1.0        1st Qu.:2.0     Class :character  
##  Mode  :character   Median :2.5        Median :3.5     Mode  :character  
##                     Mean   :2.5        Mean   :3.5                       
##                     3rd Qu.:4.0        3rd Qu.:5.0                       
##                     Max.   :5.0        Max.   :6.0                       
##                                                                          
##  Date_Start_Eggs    Date_End_Eggs       Image_Box1         Image_Box2       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     MC_Box1          MC_Box2       MC_Total_Adults     HC_Box1      
##  Min.   :  0.00   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
##  1st Qu.: 26.52   1st Qu.: 25.99   1st Qu.: 51.05   1st Qu.: 22.00  
##  Median : 58.87   Median : 54.29   Median :101.50   Median : 56.50  
##  Mean   : 75.88   Mean   : 70.57   Mean   :128.37   Mean   : 72.06  
##  3rd Qu.:106.93   3rd Qu.: 99.40   3rd Qu.:170.25   3rd Qu.:101.00  
##  Max.   :327.40   Max.   :299.00   Max.   :514.40   Max.   :288.00  
##  NA's   :142      NA's   :141      NA's   :5        NA's   :144     
##     HC_Box2       HC_Total_Adults   Nb_parents      Growth_rate    
##  Min.   :  0.00   Min.   :  0.0   Min.   :  0.00   Min.   :0.0000  
##  1st Qu.: 25.00   1st Qu.: 62.0   1st Qu.: 67.25   1st Qu.:0.4183  
##  Median : 52.00   Median :100.0   Median :100.00   Median :0.9386  
##  Mean   : 68.08   Mean   :131.8   Mean   :138.30   Mean   :1.4482  
##  3rd Qu.: 96.75   3rd Qu.:169.0   3rd Qu.:188.00   3rd Qu.:2.3587  
##  Max.   :290.00   Max.   :499.0   Max.   :499.00   Max.   :6.8571  
##  NA's   :142      NA's   :43      NA's   :122      NA's   :142     
##  First_Throw        Extinction_finalstatus Less_than_5       
##  Length:576         Length:576             Length:576        
##  Class :character   Class :character       Class :character  
##  Mode  :character   Mode  :character       Mode  :character  
##                                                              
##                                                              
##                                                              
## 
#Remove populations killed due to the pathogen
data_status  <- data.frame(data$Population,data$Extinction_finalstatus)
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)
## [1] 72 23
data<-data[data$Extinction_finalstatus!="kill",]

#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove 

#Check variables
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : chr  "Block3" "Block3" "Block3" "Block3" ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : chr  "High1" "High2" "High3" "High4" ...
##  $ Genetic_Diversity     : chr  "High" "High" "High" "High" ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: chr  "no" "no" "no" "no" ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)

#Order levels 
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)
## [1] "High"   "Medium" "Low"
data<-droplevels(data) #remove previous levels 
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
##  $ Genetic_Diversity     : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5 Nb_adults Lambda obs
## 1          no       100     NA   1
## 2          no       100     NA   2
## 3          no       100     NA   3
## 4          no       100     NA   4
## 5          no       100     NA   5
## 6          no       100     NA   6
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571            Medium                  5               6     8/25/21
## 572            Medium                  5               6     8/25/21
## 573            Medium                  5               6     8/25/21
## 574            Medium                  5               6     8/25/21
## 575            Medium                  5               6     8/25/21
## 576            Medium                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults   Lambda obs
## 571          no                     no         yes         4 0.500000 499
## 572     extinct                    yes         yes        NA       NA 500
## 573          no                     no          no        46 1.314286 501
## 574     extinct                    yes         yes        NA       NA 502
## 575          no                     no          no        56 1.400000 503
## 576          no                     no          no       133 3.243902 504
dim(data)
## [1] 504  26
#Add variable 
data$gen_square <- data$Generation_Eggs^2

1.2 Heterozygosity remain over time

#Upload He dataset
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]

#Merge dataset: add heterozygosity data to oridinal data 
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)

### Creation of 4 He: 
  # 1- He_start: He was remaining for each population when we started the experiment 
  # 2- He_lost_gen_t: He was lost only during this specfic generation during the experiment  (considering He=1 before this generation)
  # 3- He_remain: He was remaining after each generation (He_remain[Gen=6]=He_end)
  # 4- He_exp: He was lost during the entire experiment (considering He=1 at the beginning of the exp)
  # 5- He_end: He was remaining for each population when we finished the experiment (considering the real He at the beginning of the exp, He_start, not 1)


##### 
##### 1- He_start
#####
colnames(data)[28] <- "He_start"



##### 
##### 2- He_lost_gen_t
##### 
#Add He lost each generation
data$He_lost_gen_t <- 1 - (1/(2*data$Nb_adults))
#To consider extinct populations
is.na(data$He_lost_gen_t) <- sapply (data$He_lost_gen_t, is.infinite)


##### 
##### 3- He_remain
##### 
#Create new variable He remain per generation
data$He_remain <- NA
#Initiate with gen=1: He_start * He_lost_gen1
data$He_remain[data$Generation_Eggs=="1"] <- data$He_start[data$Generation_Eggs=="1"]*
  data$He_lost_gen_t[data$Generation_Eggs=="1"]

#Compute for all the other genrations, except the first one 
for(i in levels(data$Population)){
     ifelse(sum(data$Generation_Eggs[data$Population==i])!=21,print("ERROR"),print(i))
 for(j in 2:max(data$Generation_Eggs)){
   data$He_remain[data$Population==i&data$Generation_Eggs==j] <- 
     data$He_remain[data$Population==i&data$Generation_Eggs==j-1]*
     data$He_lost_gen_t[data$Population==i&data$Generation_Eggs==j]
   }
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
##### 
##### 4- He_exp
##### 
# He at the end of the experiment 
data$He_exp <- NA 

#Compute the HE at the end of the experiment 
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Eggs)!=21,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_adults))
#To consider extinct add this row: 
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
#Compute for the whole experiment 
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data$He_exp[data$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
##### 
##### 5- He_end
##### 
# Remaining He at the end of the experiment
data$He_end <- data$He_start * data$He_exp

#Save these values for Phenotyping 
data_he <- merge(x = data_he, 
                 y = data[data$Generation_Eggs==1,
                          c("Population","He_start",
                             "He_exp","He_end")], 
                 by = "Population", all.x=TRUE)

1.3 Phenotyping dataset

#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)


#Factors variables 
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <-   plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity, 
                                             levels = c("High", "Medium", "Low"))



#Add sum total 
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults  / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx

data_phenotyping_4week <- data_phenotyping[data_phenotyping$Week=="4-week",]
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]




#Merge with He dataset
data_phenotyping <- merge(x = data_phenotyping, 
              y = data_he, by = "Population", all.x=TRUE)
data_phenotyping_4week <- merge(x = data_phenotyping_4week, 
              y = data_he, by = "Population", all.x=TRUE)

#Clean if less than 30 parents 
pop_possible<-unique(data$Population[data$Generation_Eggs==5&data$Nb_adults>=30])
data_phenotyping <- data_phenotyping[data_phenotyping$Population %in% pop_possible,]

pop_possible<-unique(data$Population[data$Generation_Eggs==4&data$Nb_adults>=30])
data_phenotyping_4week <- data_phenotyping_4week[data_phenotyping_4week$Population %in% pop_possible,]

data_phenotyping_all <- rbind(data_phenotyping,data_phenotyping_4week)

1.4 Old He remain

# #Upload growth rate end of experiment
# data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
# data_he$Population <- as.factor(data_he$Population)
# data_he <- data_he[,c(1,3)]
# 
# 
# ### Additional: He remaining 
# # New vector for the lost during exp
# data_he$He_remain_exp <- NA
# 
# # He lost during exp
# for(i in levels(data$Population)){
# temp_data_pop <- subset(data, Population==i)
# ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
# temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))
# 
# #To consider extinct add this row: 
# is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
# 
# temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
# data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
# rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
# }
# 
# # Remaining He at the end of the experiment
# data_he$He_end <- data_he$He_remain * data_he$He_remain_exp
# 
# # Remove kill population 
# data_he <- data_he[!is.na(data_he$He_remain_exp),]
# data_he <- droplevels(data_he)
# 
# ## Merge dataset He
# data_old_merged <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)

1.5 Survival data

library(survival)
library(lubridate)
library(ggsurvfit)
library("gtsummary")
library("tidycmprsk")
library("condSURV")

#########


data_survival <- data[data$Generation_Parents==0,c("Population","Block","Genetic_Diversity")]
data_survival$Gen_eggs_extinct <- max(data$Generation_Eggs) 
data_survival$Survival <- "yes" 

#Format data
for (i in unique(data_survival$Population)) {
  if (data$Extinction_finalstatus[data$Population==i&data$Generation_Eggs==1] == "yes" ) {
    data_survival$Gen_eggs_extinct[data_survival$Population==i]<-
      min(data$Generation_Eggs[data$Population==i&data$First_Throw=="extinct"])
    data_survival$Survival[data_survival$Population==i] <- "extinct"
  }else{
    print(i)
  }
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low29"
## [1] "Low3"
## [1] "Low34"
## [1] "Low35"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med15"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med21"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
  str(data_survival)
## 'data.frame':    84 obs. of  5 variables:
##  $ Population       : Factor w/ 84 levels "High1","High13",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Block            : Factor w/ 3 levels "Block3","Block4",..: 1 2 2 2 2 2 2 2 1 2 ...
##  $ Genetic_Diversity: Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Gen_eggs_extinct : int  6 6 6 6 6 6 6 6 6 6 ...
##  $ Survival         : chr  "yes" "yes" "yes" "yes" ...
data_survival$Survival <- as.factor(data_survival$Survival)
#Note: the Surv() function in the {survival} package accepts by default TRUE/FALSE, where TRUE is event and FALSE is censored; 1/0 where 1 is event and 0 is censored; or 2/1 where 2 is event and 1 is censored. Please take care to ensure the event indicator is properly formatted.
data_survival$status <- ifelse(data_survival$Survival=="extinct", 1, 0)

2 PLOT

2.1 Population size

PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults,  
                                group = Population,
                                colour = Genetic_Diversity)) + 
  geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
  geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

SUM_Popsize <- Rmisc::summarySE(data, 
                        measurevar = c("Nb_adults"),
                       groupvars = c("Generation_Eggs",
                                     "Genetic_Diversity"), 
                       na.rm=TRUE)
SUM_Popsize
##    Generation_Eggs Genetic_Diversity  N Nb_adults        sd        se       ci
## 1                1              High 27 100.00000   0.00000  0.000000  0.00000
## 2                1            Medium 23 100.00000   0.00000  0.000000  0.00000
## 3                1               Low 34 100.00000   0.00000  0.000000  0.00000
## 4                2              High 27 344.29630  73.77294 14.197610 29.18360
## 5                2            Medium 23 282.04348 100.75735 21.009360 43.57075
## 6                2               Low 34 282.61765  78.53006 13.467794 27.40043
## 7                3              High 27 151.85185  80.96899 15.582489 32.03027
## 8                3            Medium 23 118.73913  88.54491 18.462891 38.28969
## 9                3               Low 34 104.61765  85.13341 14.600260 29.70445
## 10               4              High 26 114.00000  80.20224 15.728954 32.39439
## 11               4            Medium 23  62.65217  64.50483 13.450188 27.89398
## 12               4               Low 34  46.38235  39.63241  6.796903 13.82840
## 13               5              High 27 103.62963  51.56486  9.923661 20.39838
## 14               5            Medium 21  58.57143  51.11122 11.153383 23.26555
## 15               5               Low 31  51.51613  46.13377  8.285870 16.92200
## 16               6              High 27 147.66667  42.60282  8.198916 16.85311
## 17               6            Medium 19  69.73684  46.02396 10.558620 22.18284
## 18               6               Low 28  77.03571  46.14428  8.720450 17.89289
# SUM_Popsize <- Rmisc::summarySE(data[data$Extinction_finalstatus=="no",],
#                         measurevar = c("Nb_adults"),
#                        groupvars = c("Generation_Eggs",
#                                      "Genetic_Diversity"),
#                        na.rm=TRUE)
# SUM_Popsize



PLOT_Pop_size_average <- PLOT_Pop_size + 
  geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 5, position = position_dodge(0.4)) +
  geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                ymin = Nb_adults-ci, ymax = Nb_adults+ci, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 1, width=.2,  position = position_dodge(0.4)) +
  geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
           linetype = "dashed", size = 1.5, position = position_dodge(0.4)) 
PLOT_Pop_size_average

# cowplot::save_plot(file = here::here("figures","PopulationSize_average.pdf"),   PLOT_Pop_size_average, base_height = 10/cm(1),           base_width = 15/cm(1), dpi = 610)






#Prediction with models 
# Test effect

# data$gen_square <- data$Generation_Eggs^2
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity +
                      gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], 
            family = "poisson")



#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(2,6, length = 100),each=3),
                       Block = "Block4",
                       gen_square = (rep(seq(2,6, length = 100),each=3))^2,
                       Estimates = NA)

filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)


#Predcit estimate minimun values 
tapply(filldata$Estimates, filldata$Genetic_Diversity, min)
##     High      Low   Medium 
## 85.85198 44.16308 51.99642
filldata[filldata$Estimates>=44.16&filldata$Estimates<=44.17,]
##     Genetic_Diversity Generation_Eggs  Block gen_square Estimates
## 198               Low        4.626263 Block4   21.40231  44.16308
filldata[filldata$Estimates>=51.99&filldata$Estimates<=52.00,]
##     Genetic_Diversity Generation_Eggs  Block gen_square Estimates
## 230            Medium        5.070707 Block4   25.71207  51.99642
filldata[filldata$Estimates>=85.85&filldata$Estimates<=85.86,]
##     Genetic_Diversity Generation_Eggs  Block gen_square Estimates
## 181              High        4.424242 Block4   19.57392  85.85198
## Change name vector 
vector_names <- c(`Low` = "Strong bottleneck",
                    `Medium` = "Intermediate bottleneck",
                    `High` = "Diverse")


PLOT_Pop_size_predict<-ggplot2::ggplot(data[data$Extinction_finalstatus=="no",], 
                               aes(x = factor(Generation_Eggs), 
                                y = Nb_adults)) + 
  geom_point(aes(group = Population,
                 colour = Genetic_Diversity), 
             position = position_dodge(0.4), size = 0.7, alpha = 0.5) +
  geom_line(aes(group = Population,
                colour = Genetic_Diversity), 
            position = position_dodge(0.4), size = 0.05, alpha = 0.5) +
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.4) + 
  scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07", "#E7B800")) +
  ylab("Population size") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size_predict

# 
# cowplot::save_plot(file = here::here("figures","PopulationSize_predict.pdf"),   PLOT_Pop_size_predict, base_height = 10/cm(1),           base_width = 15/cm(1), dpi = 610)
# 
# 

2.2 Extinction

#Percent of extinction
length(unique(data$Population[data$Genetic_Diversity=="Low"&
                                data$Extinction_finalstatus=="yes"]))/
  length(unique(data$Population[data$Genetic_Diversity=="Low"]))
## [1] 0.2352941
length(unique(data$Population[data$Genetic_Diversity=="Medium"&
                                data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="Medium"]))
## [1] 0.2173913
length(unique(data$Population[data$Genetic_Diversity=="High"&
                                data$Extinction_finalstatus=="yes"]))/
  length(unique(data$Population[data$Genetic_Diversity=="High"]))
## [1] 0
#Create data with percent
vector_names <- c(`Low` = "Strong bottleneck",
                    `Medium` = "Intermediate bottleneck",
                    `High` = "Diverse")
f_labels = data.frame(Genetic_Diversity = c("High","Medium","Low"), 
                       label = c("Extinction = 0 %", "Extinction = 21.7 %", "Extinction = 23.5 %"))
f_labels$Genetic_Diversity <- factor(f_labels$Genetic_Diversity, 
                                     levels = c("High", "Medium", "Low"))






## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data,  aes(x = factor(Generation_Eggs), 
                                y = Nb_adults)) + 
  facet_wrap(~ Genetic_Diversity, ncol=3, labeller = ggplot2::as_labeller(vector_names)) + 
  #geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
  geom_line(aes(group = Population,
                colour = Extinction_finalstatus), 
            position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
  geom_text(x = 5, y = 420, 
            aes(label = label), 
            data = f_labels, 
            col="red", 
            size = 3, 
            vjust = 0) +
  scale_color_manual(values = c("black", "red")) +
  ylab("Population size") +
  xlab("Generation") +
  theme(legend.position = "none") +
  theme_LO
PLOT_Extinction

# 
# 
# cowplot::save_plot(here::here("figures","Extinction.pdf"),   PLOT_Extinction, base_height = 8/cm(1),
#           base_width = 20/cm(1), dpi = 610)
# #
# 



PLOT_Extinction 

2.3 Growth rate G2

tapply(data[data$Generation_Eggs==2,]$Lambda,data[data$Generation_Eggs==2,]$Genetic_Diversity,mean)
##     High   Medium      Low 
## 3.442963 2.820435 2.826176
PLOT_Growth<-ggplot2::ggplot(data[data$Generation_Eggs==2,], aes(x = Genetic_Diversity, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
  geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Genetic diversity") +
  theme_LO
PLOT_Growth

# 
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"),   PLOT_Growth, base_height = 10/cm(1),
#           base_width = 8/cm(1), dpi = 610)

2.4 Life stage proportion

SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_adults"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_adults
##   Genetic_Diversity   Week  N  p_adults         sd         se         ci
## 1              High 4-week 40 0.6315807 0.18084723 0.02859446 0.05783775
## 2              High 5-week 94 0.9646852 0.05234098 0.00539856 0.01072047
## 3            Medium 4-week 18 0.7261684 0.16721535 0.03941304 0.08315424
## 4            Medium 5-week 38 0.9594668 0.06046606 0.00980889 0.01987470
## 5               Low 4-week 27 0.6649974 0.19731065 0.03797245 0.07805350
## 6               Low 5-week 52 0.9538157 0.08241700 0.01142918 0.02294504
SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_pupae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_pupae
##   Genetic_Diversity   Week  N    p_pupae         sd          se          ci
## 1              High 4-week 40 0.31153108 0.15526911 0.024550202 0.049657471
## 2              High 5-week 94 0.01417884 0.02356334 0.002430373 0.004826238
## 3            Medium 4-week 18 0.23410139 0.15149574 0.035707888 0.075337059
## 4            Medium 5-week 38 0.01931921 0.02904804 0.004712215 0.009547854
## 5               Low 4-week 27 0.28241476 0.17824104 0.034302504 0.070509807
## 6               Low 5-week 52 0.01796037 0.02546214 0.003530964 0.007088705
SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_larvae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_larvae
##   Genetic_Diversity   Week  N   p_larvae         sd          se          ci
## 1              High 4-week 40 0.05688822 0.04383041 0.006930197 0.014017646
## 2              High 5-week 94 0.02113595 0.04275336 0.004409672 0.008756736
## 3            Medium 4-week 18 0.03973025 0.03784416 0.008919954 0.018819458
## 4            Medium 5-week 38 0.02121399 0.04478738 0.007265472 0.014721244
## 5               Low 4-week 27 0.05258789 0.08541162 0.016437473 0.033787710
## 6               Low 5-week 52 0.02822390 0.07073457 0.009809120 0.019692631
SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
                       p_pupae=SUM_Prop_pupae$p_pupae,
                       p_larvae=SUM_Prop_larvae$p_larvae)


SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop
##    Genetic_Diversity   Week  N    Stage Proportion
## 1               High 4-week 40 p_adults 0.63158069
## 2               High 5-week 94 p_adults 0.96468521
## 3             Medium 4-week 18 p_adults 0.72616835
## 4             Medium 5-week 38 p_adults 0.95946680
## 5                Low 4-week 27 p_adults 0.66499735
## 6                Low 5-week 52 p_adults 0.95381574
## 7               High 4-week 40  p_pupae 0.31153108
## 8               High 5-week 94  p_pupae 0.01417884
## 9             Medium 4-week 18  p_pupae 0.23410139
## 10            Medium 5-week 38  p_pupae 0.01931921
## 11               Low 4-week 27  p_pupae 0.28241476
## 12               Low 5-week 52  p_pupae 0.01796037
## 13              High 4-week 40 p_larvae 0.05688822
## 14              High 5-week 94 p_larvae 0.02113595
## 15            Medium 4-week 18 p_larvae 0.03973025
## 16            Medium 5-week 38 p_larvae 0.02121399
## 17               Low 4-week 27 p_larvae 0.05258789
## 18               Low 5-week 52 p_larvae 0.02822390
SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))



# New facet label names for supp variable
supp.weeks <- c("4 week", "5 week")
names(supp.weeks) <- c("4-week", "5-week")
 


PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity, 
                                                  y = Proportion, fill = Stage)) +
  facet_wrap(~ Week, ncol=2, 
             labeller = labeller(Week = supp.weeks)) +
  geom_bar(stat="identity") + 
  xlab("Genetic history") +
  labs(color="Stage of individuals") +
  ylab("Proportion of each stage") +
  scale_fill_manual(values= c("#619CFF","#F8766D","#00BA38"), 
                    breaks = c("p_adults", "p_pupae", "p_larvae"),
                    labels = c( "Adult", "Pupae","Larvae")) + 
  ggplot2::scale_x_discrete(breaks=c("High", "Medium", "Low"), 
                            labels=c("Diverse", 
                                     "Intermediate\nbottleneck", 
                                      "Strong \nbottleneck")) + 
  theme_LO + 
  theme(strip.background = element_rect(color="grey30", 
                                        fill="white", size=0.5, linetype="solid"), 
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))
PLOT_prop

# 
# cowplot::save_plot(here::here("figures","Life_stage_proportion.pdf"),   PLOT_prop,
#           base_height = 10/cm(1), base_width = 18/cm(1), dpi = 610)

2.5 Growth rate and He

# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_phenotyping,
                            measurevar = c("Lambda"),
                            groupvars = c("Genetic_Diversity",
                                          "He_end",
                                          "Population"), 
                            na.rm=TRUE)


PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_end,
                                                      y = Lambda, 
                                                      colour = Genetic_Diversity, 
                                                      size = N)) + 
    geom_point(alpha = 0.8) +
    ylab("Growth rate") +
    xlab("Expected heterozygosity") +
    scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07","#E7B800")) +
    labs(size="Replicates") + 
    theme_LO 
PLOT_He_Final

# 
# 
# cowplot::save_plot(here::here("figures","Heterozygosity_Growthrate.pdf"),
#                    PLOT_He_Final,
#           base_height = 12/cm(1), base_width = 16/cm(1), dpi = 610)
# 
# 



# # ALL WITHOUT PSEUDO
# SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
#                             measurevar = c("Lambda"),
#                        groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_end",
#                                      "Population"), 
#                        na.rm=TRUE)
# 
# SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_end)
# 
# SUM_POP_He_ALL <- SUM_POP_He_ALL[!is.na(SUM_POP_He_ALL$He_end),]
# SUM_POP_He_ALL <- SUM_POP_He_ALL[SUM_POP_He_ALL$HE!="-Inf",]
# 
# 
# facet_names <- c('Initial'="Before rescue experiment",'Final'="After rescue experiment")
# 
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
#                          colour = Genetic_Diversity, 
#                          shape = Phenotyping)) +
#   #facet_wrap(~ Phenotyping, ncol=1) +
#   geom_point(size = 3, alpha = 0.8) +
#   scale_shape_manual(values = c(1, 19)) + 
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO 
# PLOT_He_ALL
# 
# 
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
#                          colour = Genetic_Diversity, 
#                          shape = Phenotyping)) +
#   facet_wrap(~ Phenotyping, ncol=1) +
#   geom_point(size = 3, alpha = 0.8) +
#   scale_shape_manual(values = c(1, 19)) + 
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO
# PLOT_He_ALL
# 
# 
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
#                          colour = Genetic_Diversity)) +
#   facet_wrap(~ Phenotyping, ncol=1) +
#   geom_point(size = 3, alpha = 0.8) +
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO
# PLOT_He_ALL

2.6 Lifestage and He

SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_adults"),
                       groupvars = c("Genetic_Diversity","Week", "Population", "He_end"), 
                       na.rm = TRUE)
colnames(SUM_Prop_adults)[6]<-"prop" 
SUM_Prop_adults$lifestage <- "adults" 


SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_pupae"),
                       groupvars = c("Genetic_Diversity","Week", "Population", "He_end"), 
                       na.rm = TRUE)
colnames(SUM_Prop_pupae)[6]<-"prop" 
SUM_Prop_pupae$lifestage <- "pupae" 

SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_larvae"),
                       groupvars = c("Genetic_Diversity","Week", 
                                     "Population", "He_end"), 
                       na.rm = TRUE)
colnames(SUM_Prop_larvae)[6]<-"prop" 
SUM_Prop_larvae$lifestage <- "larvae" 


SUM_Prop_POP <- rbind(SUM_Prop_adults, SUM_Prop_pupae, SUM_Prop_larvae)



PLOT_Proportion_lifestage <- ggplot2::ggplot(SUM_Prop_POP, aes(x = He_end,
                                                      y = prop, 
                                                      colour = Genetic_Diversity, 
                                                      size = N)) + 
    facet_grid( Week ~ factor(lifestage, levels=c('larvae','pupae','adults'))) +
    geom_point(alpha = 0.8) +
    ylab("Proportion") +
    xlab("Expected heterozygosity") +
    scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07", "#E7B800")) + 
    labs(size="Replicates") + 
    theme_LO 
PLOT_Proportion_lifestage

 # 
 # cowplot::save_plot(here::here("figures","Life_stage_proportion_He.pdf"),
 #            PLOT_Proportion_lifestage,
 #           base_height = 14/cm(1), base_width = 36/cm(1), dpi = 610)
 # 

# 

2.7 He over time

#Compute for all the other genrations, except the first one 
data$ID_extinction <- "extant"

for(i in unique(data$Population[data$Extinction_finalstatus=="yes"])){
    gen_all <- data$Generation_Eggs[data$Population==i&!is.na(data$He_remain)]
    maxgen <- max(gen_all)
    data$ID_extinction[data$Population==i&data$Generation_Eggs==maxgen] <- "willextinct" 
}

data[data$ID_extinction=="willextinct",]
##     Population  Block         Old_Label                  Label
## 205      Low16 Block4             B4-O1 7/14 BE B4 | G5 Low 16
## 273      Low27 Block5             B5-B1 6/16 BE B5 | G4 Low 27
## 278      Low28 Block5             B5_E1 5/12 BE B5 | G3 Low 28
## 299      Low30 Block5             B5-D1 6/16 BE B5 | G4 Low 30
## 303      Low31 Block5         B(2)5-P1  6/16 BE B5 | G4 Low 31
## 310      Low32 Block5         B(2)5_Q1  5/12 BE B5 | G3 Low 32
## 317      Low33 Block5         B(2)5-T1  7/21 BE B5 | G5 Low 33
## 336      Low36 Block5         B(2)5-X1  6/16 BE B5 | G4 Low 36
## 402      Med14 Block4   B4-Backup Fam L  6/9 BE B4 | G4 Med 14
## 409      Med17 Block5   B5_Backup_Fam_F 5/12 BE B5 | G3 Med 17
## 437      Med20 Block5 B5 - Backup Fam I 6/16 BE B5 | G4 Med 20
## 449      Med22 Block5   B5_Backup_Fam_B 5/12 BE B5 | G3 Med 22
## 479       Med5 Block3 B3 - Backup Fam E   7/7 BE B3 | G5 Med 5
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 205               Low                  4               5     7/14/21
## 273               Low                  3               4     6/16/21
## 278               Low                  2               3     5/12/21
## 299               Low                  3               4     6/16/21
## 303               Low                  3               4     6/16/21
## 310               Low                  2               3     5/12/21
## 317               Low                  4               5     7/21/21
## 336               Low                  3               4     6/16/21
## 402            Medium                  3               4      6/9/21
## 409            Medium                  2               3     5/12/21
## 437            Medium                  3               4     6/16/21
## 449            Medium                  2               3     5/12/21
## 479            Medium                  4               5      7/7/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2   MC_Box1   MC_Box2
## 205         7/14/21       7/15/21       <NA>   DSC_0994  0.000000  2.000000
## 273         6/16/21       6/17/21   DSC_0526   DSC_0527  5.090878  0.000000
## 278         5/12/21       5/13/21   DSC_0918   DSC_0919  3.085768  3.638847
## 299         6/16/21       6/17/21   DSC_0559       <NA>  1.383920  0.000000
## 303         6/16/21       6/17/21       <NA>   DSC_0541  0.000000  1.000000
## 310         5/12/21       5/13/21   DSC_0953   DSC_0954 88.254170 70.952124
## 317         7/21/21       7/22/21   DSC_0117   DSC_0118  1.000000  1.000000
## 336         6/16/21       6/17/21   DSC_0540       <NA>  1.000000  0.000000
## 402          6/9/21       6/10/21   DSC_0376   DSC_0377 10.974990  1.000000
## 409         5/12/21       5/13/21   DSC_0916   DSC_0917  1.543122  1.084327
## 437         6/16/21       6/17/21   DSC_0542       <NA>  1.000000  0.000000
## 449         5/12/21       5/13/21   DSC_0922   DSC_0923  9.168447  6.251069
## 479          7/7/21        7/8/21       <NA>   DSC_0830        NA        NA
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 205             2.0       0       2               2          5 0.400000000
## 273             5.1       3       0               3         19 0.157894737
## 278             6.7       1       1               2        374 0.005347594
## 299             1.4       1       0               1        146 0.006849315
## 303             1.0       1       1               2        272 0.007352941
## 310           159.2      71      56             127        276 0.460144928
## 317             2.0       1       1               2          3 0.666666667
## 336             1.0       1       0               1         22 0.045454545
## 402            12.0       7       1               8        138 0.057971014
## 409             2.6       3       2               5        416 0.012019231
## 437             1.0       1       0               1         20 0.050000000
## 449            15.4      10       6              16        200 0.080000000
## 479             0.0      NA      NA               1         11 0.090909091
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults      Lambda obs
## 205          no                    yes         yes         2 0.400000000 378
## 273          no                    yes         yes         3 0.157894737 319
## 278          no                    yes         yes         2 0.005347594 236
## 299          no                    yes         yes         1 0.006849315 322
## 303          no                    yes         yes         2 0.007352941 323
## 310          no                    yes         yes       127 0.460144928 240
## 317          no                    yes         yes         2 0.666666667 409
## 336          no                    yes         yes         1 0.045454545 328
## 402          no                    yes         yes         8 0.057971014 308
## 409          no                    yes         yes         5 0.012019231 245
## 437          no                    yes         yes         1 0.050000000 332
## 449          no                    yes         yes        16 0.080000000 250
## 479          no                    yes         yes         1 0.090909091 359
##     gen_square  He_start He_lost_gen_t He_remain    He_exp    He_end
## 205         25 0.5544299     0.7500000 0.3575672 0.6449276 0.3575672
## 273         16 0.5367998     0.8333333 0.4324691 0.8056432 0.4324691
## 278          9 0.5386585     0.7500000 0.4014365 0.7452523 0.4014365
## 299         16 0.5540444     0.5000000 0.2742819 0.4950540 0.2742819
## 303         16 0.5532775     0.7500000 0.4116634 0.7440450 0.4116634
## 310          9 0.5528880     0.9960630 0.5469651 0.9892872 0.5469651
## 317         25 0.5523333     0.7500000 0.3359893 0.6083089 0.3359893
## 336         16 0.5427371     0.5000000 0.2635269 0.4855518 0.2635269
## 402         16 0.6802331     0.9375000 0.6315974 0.9285014 0.6315974
## 409          9 0.6825509     0.9000000 0.6104897 0.8944237 0.6104897
## 437         16 0.6819650     0.5000000 0.3302071 0.4841994 0.3302071
## 449          9 0.6809168     0.9687500 0.6546991 0.9614965 0.6546991
## 479         25 0.6807065     0.5000000 0.3209456 0.4714889 0.3209456
##     ID_extinction
## 205   willextinct
## 273   willextinct
## 278   willextinct
## 299   willextinct
## 303   willextinct
## 310   willextinct
## 317   willextinct
## 336   willextinct
## 402   willextinct
## 409   willextinct
## 437   willextinct
## 449   willextinct
## 479   willextinct
dim(data[data$ID_extinction=="willextinct",])
## [1] 13 33
data[data$ID_extinction=="willextinct","Population"]
##  [1] Low16 Low27 Low28 Low30 Low31 Low32 Low33 Low36 Med14 Med17 Med20 Med22
## [13] Med5 
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9
data$Population[data$Genetic_Diversity=="Medium"&data$Extinction_finalstatus=="yes"]
##  [1] Med14 Med14 Med14 Med14 Med14 Med14 Med17 Med17 Med17 Med17 Med17 Med17
## [13] Med20 Med20 Med20 Med20 Med20 Med20 Med22 Med22 Med22 Med22 Med22 Med22
## [25] Med5  Med5  Med5  Med5  Med5  Med5 
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9
PLOT_He_extinction <-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = He_remain, 
                                group = Population,
                                colour = Extinction_finalstatus, 
                                shape = ID_extinction)) + 
  geom_line(position = position_dodge(0.5), 
            size = 0.25, alpha = 0.6) +
  geom_point(position = position_dodge(0.5), size = 3,  stroke = 0.6,
             fill= "white", alpha = 0.8) +
  ylab("Expected heterozygozity") +
  scale_color_manual(values = c("black", "red")) +
  scale_shape_manual(values = c(21, 13), guide=FALSE) +
  labs(color="Extinct populations") +
  xlab("Generation") +
  theme_LO
PLOT_He_extinction

# # 
# cowplot::save_plot(here::here("figures","Heterozygosity_over_time.pdf"),
#                    PLOT_He_extinction,
#           base_height = 12/cm(1), base_width = 14/cm(1), dpi = 610)
# 

2.8 Survival vs He

s2 <- survfit(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
summary(s2)
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, 
##     data = data_survival)
## 
##                 Genetic_Diversity=High 
##      time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 
##                 Genetic_Diversity=Medium 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     23       2    0.913  0.0588        0.805        1.000
##     5     21       2    0.826  0.0790        0.685        0.996
##     6     19       1    0.783  0.0860        0.631        0.971
## 
##                 Genetic_Diversity=Low 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     34       2    0.941  0.0404        0.865        1.000
##     5     32       4    0.824  0.0654        0.705        0.962
##     6     28       2    0.765  0.0727        0.635        0.921
# Create a summary dataset 
SUM_survprob <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity), each =6), 
                       Generation_Eggs = rep(seq(1:6), n=3),
                       SurvProb = 1, 
                       se = 0)


SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Medium"&
                        SUM_survprob$Generation_Eggs==4]<-summary(s2)$surv[1]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Medium"&
                        SUM_survprob$Generation_Eggs==5]<-summary(s2)$surv[2]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Medium"&
                        SUM_survprob$Generation_Eggs==6]<-summary(s2)$surv[3]

SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Low"&
                        SUM_survprob$Generation_Eggs==4]<-summary(s2)$surv[4]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Low"&
                        SUM_survprob$Generation_Eggs==5]<-summary(s2)$surv[5]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Low"&
                        SUM_survprob$Generation_Eggs==6]<-summary(s2)$surv[6]
  
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Medium"&
                        SUM_survprob$Generation_Eggs==4]<-summary(s2)$std.err[1]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Medium"&
                        SUM_survprob$Generation_Eggs==5]<-summary(s2)$std.err[2]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Medium"&
                        SUM_survprob$Generation_Eggs==6]<-summary(s2)$std.err[3]

SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Low"&
                        SUM_survprob$Generation_Eggs==4]<-summary(s2)$std.err[4]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Low"&
                        SUM_survprob$Generation_Eggs==5]<-summary(s2)$std.err[5]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Low"&
                        SUM_survprob$Generation_Eggs==6]<-summary(s2)$std.err[6]
  

data_survprob_he <- merge(x = data, y = SUM_survprob, by = c("Genetic_Diversity", "Generation_Eggs"), all.x=TRUE)
head(data_survprob_he)
##   Genetic_Diversity Generation_Eggs Population  Block       Old_Label
## 1              High               1      High1 Block3 B3_All_Red_Mix 
## 2              High               1     High26 Block5 B5_All_Red_Mix 
## 3              High               1     High14 Block4 B4_All_Red_Mix 
## 4              High               1     High33 Block5 B5_All_Red_Mix 
## 5              High               1      High4 Block3 B3_All_Red_Mix 
## 6              High               1     High19 Block4 B4_All_Red_Mix 
##                     Label Generation_Parents Date_Census Date_Start_Eggs
## 1  BE B3 2/17 | G1  High1                  0     2/17/21         2/17/21
## 2  BE B5 3/3 | G1  High26                  0      3/3/21          3/3/21
## 3 BE B4 2/24 | G1  High14                  0     2/24/21         2/24/21
## 4  BE B5 3/3 | G1  High33                  0      3/3/21          3/3/21
## 5  BE B3 2/17 | G1  High4                  0     2/17/21         2/17/21
## 6 BE B4 2/24 | G1  High19                  0     2/24/21         2/24/21
##   Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1
## 1       2/18/21   DSC_0155       <NA>      NA      NA           101.9      NA
## 2        3/4/21   DSC_0675       <NA>      NA      NA           101.5      NA
## 3       2/25/21   DSC_0416       <NA>      NA      NA           107.5      NA
## 4        3/4/21   DSC_0682       <NA>      NA      NA            82.5      NA
## 5       2/18/21   DSC_0158       <NA>      NA      NA           101.9      NA
## 6       2/25/21   DSC_0422       <NA>      NA      NA           101.8      NA
##   HC_Box2 HC_Total_Adults Nb_parents Growth_rate First_Throw
## 1      NA             100         NA          NA          no
## 2      NA             100         NA          NA          no
## 3      NA             100         NA          NA          no
## 4      NA             100         NA          NA          no
## 5      NA             100         NA          NA          no
## 6      NA             100         NA          NA          no
##   Extinction_finalstatus Less_than_5 Nb_adults Lambda obs gen_square  He_start
## 1                     no          no       100     NA   1          1 0.9986008
## 2                     no          no       100     NA  59          1 0.9986008
## 3                     no          no       100     NA  28          1 0.9986008
## 4                     no          no       100     NA  63          1 0.9986008
## 5                     no          no       100     NA   4          1 0.9986008
## 6                     no          no       100     NA  33          1 0.9986008
##   He_lost_gen_t He_remain    He_exp    He_end ID_extinction SurvProb se
## 1         0.995 0.9936078 0.9679591 0.9666047        extant        1  0
## 2         0.995 0.9936078 0.9220546 0.9207645        extant        1  0
## 3         0.995 0.9936078 0.9766045 0.9752381        extant        1  0
## 4         0.995 0.9936078 0.9358965 0.9345870        extant        1  0
## 5         0.995 0.9936078 0.9812175 0.9798447        extant        1  0
## 6         0.995 0.9936078 0.9804728 0.9791010        extant        1  0
Plot_survival <- ggplot2::ggplot(data_survprob_he, aes(x = He_remain,  y = SurvProb ,  colour = Genetic_Diversity)) + 
    geom_point(alpha = 0.8) +
    ylab("Survival probability") +
    xlab("Expected heterozygosity") +
    scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07", "#E7B800")) + 
    labs(size="Replicates") + 
    theme_LO 

# 
# 
# cowplot::save_plot(here::here("figures","Survival_function_He.pdf"),
#                    Plot_survival,
#           base_height = 12/cm(1), base_width = 14/cm(1), dpi = 610)
# 

3 ANALYSIS

3.1 Probability of extinction

############ 
###### Clean dataset
############ 
#Prepare clean dataset
data_proba_extinction <- data[data$Generation_Eggs==6,
                              c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)



############ 
###### Analysis
############ 
# #Analysis
# mod0 <- glm(y ~ Genetic_Diversity + Block,  
#       data = data_proba_extinction, family = binomial)
# 
# mod0 <- glm(y ~ Genetic_Diversity-1 + Block,  
#       data = data_proba_extinction, family = binomial)
# 
# 
# #But issue of convergence 
# plogis(confint(mod0))


# ############ 
# ###### Analysis
# ############ 
# 
# 
# summary(mod0)
# 
# 
# 
# #Test convergence 
# mod1 <- glm(y ~ Genetic_Diversity + Block ,  
#       data = data_proba_extinction[data_proba_extinction$Genetic_Diversity!="High",], 
#       family = binomial)
# 
# #Test genetic diversity effect
# mod0 <- glm(y ~ Block ,  
#       data = data_proba_extinction[data_proba_extinction$Genetic_Diversity!="High",], family = binomial)
# anova(mod0, mod1, test = "Chisq")
# #We remove genetic diversity 
# 
# # Perform the lrtest 
# (A <- logLik(mod1))
# (B <- logLik(mod0))
# #Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
# (teststat <- -2 * (as.numeric(A)-as.numeric(B)))
# #df = 5 - 3 = 2
# (p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))
# 
# lmtest::lrtest(mod1, mod0)




############ 
###### CONFIDENCE INTERVAL 
############
# mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,  
#       data = data_proba_extinction[data_proba_extinction$Genetic_Diversity!="High",], family = binomial)
# summary(mod0)
# 
# #Extract probability of extinction 
# #(p_high <- plogis(-21.3144))
# (p_med <- plogis(-3.0142))
# (p_low <- plogis(-2.8082))
# 
# 
# # Confidence interval
# table_raw <- confint(mod0)
# table_ci <- plogis(table_raw)
# table_ci
# 
# ############ 
# ###### Predicted value
# ############ 
# 
# 
# data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
#                                    Block = "Block4")
# 
# data_predict_extinct$predict <- predict(mod0, type = "response",
#         re.form = NA,
#         newdata = data_predict_extinct)
#         
# data_predict_extinct
# p_high
# p_med
# p_low
# 
# 
# ############ 
# ###### Posthoc
# ############ 
# emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
# 






############## 
######## 1- REAL DATA: issue with high no extinction 
############## 

mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,  
      data = data_proba_extinction, family = binomial)
summary(mod0)
## 
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial, 
##     data = data_proba_extinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.26011  -0.44258  -0.30957  -0.00005   2.47474  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## Genetic_DiversityHigh    -21.3144  1855.7435  -0.011  0.99084   
## Genetic_DiversityMedium   -3.0142     1.1293  -2.669  0.00760 **
## Genetic_DiversityLow      -2.8082     1.0659  -2.635  0.00843 **
## BlockBlock4                0.7402     1.2714   0.582  0.56046   
## BlockBlock5                3.0006     1.1265   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.449  on 84  degrees of freedom
## Residual deviance:  46.833  on 79  degrees of freedom
## AIC: 56.833
## 
## Number of Fisher Scoring iterations: 18
(p_high <- plogis(-21.3144))
## [1] 5.536989e-10
(p_med <- plogis(-3.0142))
## [1] 0.04678847
(p_low <- plogis(-2.8082))
## [1] 0.05688267
# Confidence interval: high [0,1]
plogis(confint(mod0))
##                                 2.5 %    97.5 %
## Genetic_DiversityHigh   1.689103e-236 1.0000000
## Genetic_DiversityMedium  2.415647e-03 0.2298311
## Genetic_DiversityLow     3.204498e-03 0.2431697
## BlockBlock4              1.553712e-01 0.9794348
## BlockBlock5              7.579725e-01 0.9975064
## ISSUE: DOESNT CONVERGE



mod0 <- glm(y ~ Genetic_Diversity + Block ,  
      data = data_proba_extinction, family = binomial)
mod1 <- glm(y ~ Block ,  
      data = data_proba_extinction, family = binomial)
anova(mod1,mod0,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1        81     58.903                        
## 2        79     46.833  2   12.069 0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############## 
######## 2- Augmented Data: As in Scucz et al. 2017 - Ecology
############## 
data_proba_extinction_augmented <- rbind(data_proba_extinction,
                                         # c("Block3","HighX","High","no",0), 
                                         # c("Block3","LowX","Low","no",0), 
                                         # c("Block3","MedX","Medium","no",0), 
                                         #                                                                                      c("Block3","LowY","Low","yes",1), 
#                                         c("Block3","MedXY","Medium","yes",1),
                                         c("Block3","HighY","High","yes",1))

data_proba_extinction_augmented$y <- as.numeric(data_proba_extinction_augmented$y)

mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,  
      data = data_proba_extinction_augmented, family = binomial)
summary(mod0)
## 
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial, 
##     data = data_proba_extinction_augmented)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1885  -0.4686  -0.4349  -0.1538   2.9692  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## Genetic_DiversityHigh   -4.39594    1.28253  -3.428 0.000609 ***
## Genetic_DiversityMedium -2.31099    0.86770  -2.663 0.007737 ** 
## Genetic_DiversityLow    -2.11810    0.79597  -2.661 0.007790 ** 
## BlockBlock4             -0.03575    1.05172  -0.034 0.972886    
## BlockBlock5              2.14410    0.86309   2.484 0.012984 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 117.835  on 85  degrees of freedom
## Residual deviance:  58.406  on 80  degrees of freedom
## AIC: 68.406
## 
## Number of Fisher Scoring iterations: 6
#Real data
(p_high <- plogis(-21.3144))
## [1] 5.536989e-10
(p_med <- plogis(-3.0142))
## [1] 0.04678847
(p_low <- plogis(-2.8082))
## [1] 0.05688267
#Augmented data
(p_high <- plogis(-3.9197))
## [1] 0.01946081
(p_med <- plogis(-1.6788))
## [1] 0.1572544
(p_low <- plogis(-1.5627))
## [1] 0.1732596
#Augmented data with only extinction
(p_high <- plogis(-4.39594))
## [1] 0.01217718
(p_med <- plogis(-2.31099))
## [1] 0.09021685
(p_low <- plogis(-2.11810))
## [1] 0.10735
# Confidence interval: high [0,1]
plogis(confint(mod0))
##                                2.5 %     97.5 %
## Genetic_DiversityHigh   0.0004918074 0.09139237
## Genetic_DiversityMedium 0.0127711273 0.30480717
## Genetic_DiversityLow    0.0174349180 0.31771819
## BlockBlock4             0.0963225771 0.89727408
## BlockBlock5             0.6478312906 0.98431182
mod0 <- glm(y ~ Genetic_Diversity + Block ,  
      data = data_proba_extinction_augmented, family = binomial)
mod1 <- glm(y ~ Block ,  
      data = data_proba_extinction_augmented, family = binomial)
anova(mod1,mod0,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1        82     64.685                       
## 2        80     58.406  2   6.2787  0.04331 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############## 
######## 3- PenaliZed  Logistic Regression data 
############## 


fit<-logistf::logistf(y ~ Genetic_Diversity-1 + Block, 
                      data=data_proba_extinction)
summary(fit)
## logistf::logistf(formula = y ~ Genetic_Diversity - 1 + Block, 
##     data = data_proba_extinction)
## 
## Model fitted by Penalized ML
## Coefficients:
##                               coef  se(coef)  lower 0.95 upper 0.95      Chisq
## Genetic_DiversityHigh   -5.5403989 1.6381117 -10.6142028 -2.9271591 31.0874311
## Genetic_DiversityMedium -2.5593086 0.9220804  -4.9074011 -0.9845788 12.1477934
## Genetic_DiversityLow    -2.3929983 0.8625221  -4.6407993 -0.9400804 12.5897375
## BlockBlock4              0.5511409 1.0469817  -1.5552881  3.0070956  0.2677213
## BlockBlock5              2.5551098 0.9273661   0.9131677  4.8771647 10.2850062
##                                    p method
## Genetic_DiversityHigh   2.466634e-08      2
## Genetic_DiversityMedium 4.914598e-04      2
## Genetic_DiversityLow    3.878706e-04      2
## BlockBlock4             6.048644e-01      2
## BlockBlock5             1.341156e-03      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=29.25035 on 5 df, p=2.070724e-05, n=84
## Wald test = 23.62508 on 5 df, p = 0.0002562504
nobs(fit)
## [1] 84
# Test genetic diversity effect
mod0<-logistf::logistf(y ~ Genetic_Diversity + Block, 
                      data=data_proba_extinction)
drop1(mod0)
##                       ChiSq df     P-value
## Genetic_Diversity  9.394033  2 0.009122452
## Block             12.698693  2 0.001747889
extractAIC(mod0)
##              null 
##  4.00000 29.29027
mod0<-logistf::logistf(y ~ Genetic_Diversity + Block, 
                      data=data_proba_extinction)
mod1<-logistf::logistf(y ~ Block, 
                      data=data_proba_extinction)
anova(mod1,mod0)
## Comparison of logistf models:
##                         Formula ChiSquared 
## 1 y ~ Genetic_Diversity + Block   21.29027 
## 2                     y ~ Block   11.89623 
## 
## Method:  nested 
## Chi-Squared:  9.394033   df= 2   P= 0.009122452
#Confirmation with LRT





### Estimates and CI 
fit<-logistf::logistf(y ~ Genetic_Diversity-1 + Block, 
                      data=data_proba_extinction)
sum_p_extinct <- summary(fit)
## logistf::logistf(formula = y ~ Genetic_Diversity - 1 + Block, 
##     data = data_proba_extinction)
## 
## Model fitted by Penalized ML
## Coefficients:
##                               coef  se(coef)  lower 0.95 upper 0.95      Chisq
## Genetic_DiversityHigh   -5.5403989 1.6381117 -10.6142028 -2.9271591 31.0874311
## Genetic_DiversityMedium -2.5593086 0.9220804  -4.9074011 -0.9845788 12.1477934
## Genetic_DiversityLow    -2.3929983 0.8625221  -4.6407993 -0.9400804 12.5897375
## BlockBlock4              0.5511409 1.0469817  -1.5552881  3.0070956  0.2677213
## BlockBlock5              2.5551098 0.9273661   0.9131677  4.8771647 10.2850062
##                                    p method
## Genetic_DiversityHigh   2.466634e-08      2
## Genetic_DiversityMedium 4.914598e-04      2
## Genetic_DiversityLow    3.878706e-04      2
## BlockBlock4             6.048644e-01      2
## BlockBlock5             1.341156e-03      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=29.25035 on 5 df, p=2.070724e-05, n=84
## Wald test = 23.62508 on 5 df, p = 0.0002562504
# Confidence interval:
plogis(sum_p_extinct$coefficients)
##   Genetic_DiversityHigh Genetic_DiversityMedium    Genetic_DiversityLow 
##             0.003909616             0.071803609             0.083708170 
##             BlockBlock4             BlockBlock5 
##             0.634400249             0.927916044
plogis(sum_p_extinct$ci.lower)
##   Genetic_DiversityHigh Genetic_DiversityMedium    Genetic_DiversityLow 
##            2.456403e-05            7.337438e-03            9.557749e-03 
##             BlockBlock4             BlockBlock5 
##            1.743238e-01            7.136479e-01
plogis(sum_p_extinct$ci.upper)
##   Genetic_DiversityHigh Genetic_DiversityMedium    Genetic_DiversityLow 
##              0.05082721              0.27198420              0.28088410 
##             BlockBlock4             BlockBlock5 
##              0.95289365              0.99243902
# Probability 
(p_high <- plogis(-5.5403989))
## [1] 0.003909616
(p_med <- plogis(-2.5593086))
## [1] 0.07180361
(p_low <- plogis(-2.3929983))
## [1] 0.08370817
 plogis(confint(fit))
##                            Lower 95%  Upper 95%
## Genetic_DiversityHigh   2.456403e-05 0.05082721
## Genetic_DiversityMedium 7.337438e-03 0.27198420
## Genetic_DiversityLow    9.557749e-03 0.28088410
## BlockBlock4             1.743238e-01 0.95289365
## BlockBlock5             7.136479e-01 0.99243902

3.2 Time to extinction

data_timeextinction<-data[data$First_Throw=="extinct",]
data_timeextinction <- data_timeextinction[complete.cases(data_timeextinction$Nb_adults), ]
data_timeextinction
##     Population  Block         Old_Label                  Label
## 207      Low16 Block4             B4-O1 8/18 BE B4 | G6 Low 16
## 274      Low27 Block5             B5-B1 7/21 BE B5 | G5 Low 27
## 277      Low28 Block5             B5-E1 6/16 BE B5 | G4 Low 28
## 296      Low30 Block5             B5-D1 7/21 BE B5 | G5 Low 30
## 301      Low31 Block5         B(2)5-P1  7/21 BE B5 | G5 Low 31
## 309      Low32 Block5         B(2)5-Q1  6/16 BE B5 | G4 Low 32
## 318      Low33 Block5         B(2)5-T1  8/25 BE B5 | G6 Low 33
## 335      Low36 Block5         B(2)5-X1  7/21 BE B5 | G5 Low 36
## 397      Med14 Block4   B4-Backup Fam L 7/14 BE B4 | G5 Med 14
## 410      Med17 Block5 B5 - Backup Fam F 6/16 BE B5 | G4 Med 17
## 438      Med20 Block5 B5 - Backup Fam I 7/21 BE B5 | G5 Med 20
## 448      Med22 Block5   B5-Backup Fam B 6/16 BE B5 | G4 Med 22
## 476       Med5 Block3 B3 - Backup Fam E  8/11 BE B3 | G6 Med 5
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 207               Low                  5               6     8/18/21
## 274               Low                  4               5     7/21/21
## 277               Low                  3               4     6/16/21
## 296               Low                  4               5     7/21/21
## 301               Low                  4               5     7/21/21
## 309               Low                  3               4     6/16/21
## 318               Low                  5               6     8/25/21
## 335               Low                  4               5     7/21/21
## 397            Medium                  4               5     7/14/21
## 410            Medium                  3               4     6/16/21
## 438            Medium                  4               5     7/21/21
## 448            Medium                  3               4     6/16/21
## 476            Medium                  5               6     8/11/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 207         8/18/21       8/19/21       <NA>       <NA>      NA      NA
## 274         7/21/21       7/22/21       <NA>       <NA>       0       0
## 277         6/16/21       6/17/21       <NA>       <NA>       0       0
## 296         7/21/21       7/22/21       <NA>       <NA>       0       0
## 301         7/21/21       7/22/21       <NA>       <NA>       0       0
## 309         6/16/21       6/17/21       <NA>       <NA>       0       0
## 318         8/25/21       8/26/21       <NA>       <NA>      NA      NA
## 335         7/21/21       7/22/21       <NA>       <NA>       0       0
## 397         7/14/21       7/15/21       <NA>       <NA>       0       0
## 410         6/16/21       6/17/21       <NA>       <NA>       0       0
## 438         7/21/21       7/22/21       <NA>       <NA>       0       0
## 448         6/16/21       6/17/21       <NA>       <NA>       0       0
## 476         8/11/21       8/12/21       <NA>       <NA>       0       0
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 207               0       0       0               0          2           0
## 274               0       0       0               0          3           0
## 277               0       0       0               0          2           0
## 296               0       0       0               0          1           0
## 301               0       0       0               0          2           0
## 309               0       0       0               0        127           0
## 318               0      NA      NA               0          2           0
## 335               0       0       0               0          1           0
## 397               0      NA      NA               0          8          NA
## 410               0       0       0               0          5           0
## 438               0       0       0               0          1           0
## 448               0       0       0               0         16           0
## 476               0      NA      NA               0          1           0
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 207     extinct                    yes         yes         0      0 462
## 274     extinct                    yes         yes         0      0 403
## 277     extinct                    yes         yes         0      0 320
## 296     extinct                    yes         yes         0      0 406
## 301     extinct                    yes         yes         0      0 407
## 309     extinct                    yes         yes         0      0 324
## 318     extinct                    yes         yes         0      0 493
## 335     extinct                    yes         yes         0      0 412
## 397     extinct                    yes         yes         0      0 392
## 410     extinct                    yes         yes         0      0 329
## 438     extinct                    yes         yes         0      0 416
## 448     extinct                    yes         yes         0      0 334
## 476     extinct                    yes         yes         0      0 443
##     gen_square  He_start He_lost_gen_t He_remain    He_exp    He_end
## 207         36 0.5544299            NA        NA 0.6449276 0.3575672
## 274         25 0.5367998            NA        NA 0.8056432 0.4324691
## 277         16 0.5386585            NA        NA 0.7452523 0.4014365
## 296         25 0.5540444            NA        NA 0.4950540 0.2742819
## 301         25 0.5532775            NA        NA 0.7440450 0.4116634
## 309         16 0.5528880            NA        NA 0.9892872 0.5469651
## 318         36 0.5523333            NA        NA 0.6083089 0.3359893
## 335         25 0.5427371            NA        NA 0.4855518 0.2635269
## 397         25 0.6802331            NA        NA 0.9285014 0.6315974
## 410         16 0.6825509            NA        NA 0.8944237 0.6104897
## 438         25 0.6819650            NA        NA 0.4841994 0.3302071
## 448         16 0.6809168            NA        NA 0.9614965 0.6546991
## 476         36 0.6807065            NA        NA 0.4714889 0.3209456
##     ID_extinction
## 207        extant
## 274        extant
## 277        extant
## 296        extant
## 301        extant
## 309        extant
## 318        extant
## 335        extant
## 397        extant
## 410        extant
## 438        extant
## 448        extant
## 476        extant
tapply(data_timeextinction$Generation_Eggs,data_timeextinction$Genetic_Diversity,mean)
##   High Medium    Low 
##     NA    4.8    5.0
#Test time of extinction
mod1 <- glm(Generation_Eggs ~ Genetic_Diversity + Block,
                    data = data_timeextinction, family = "poisson")
mod0 <- glm(Generation_Eggs ~  Block,
                    data = data_timeextinction, family = "poisson")

anova(mod0, mod1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        10    0.95459                     
## 2         9    0.74888  1  0.20571   0.6502
#We keep genetic diversity 

# Perform the lrtest 
(A <- logLik(mod1))
## 'log Lik.' -22.83384 (df=4)
(B <- logLik(mod0))
## 'log Lik.' -22.9367 (df=3)
#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))
## [1] -0.2057062
#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))
## [1] 1
lmtest::lrtest(mod0,mod1)
## Likelihood ratio test
## 
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ Genetic_Diversity + Block
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   3 -22.937                     
## 2   4 -22.834  1 0.2057     0.6502
mod1 <- glm(Generation_Eggs ~ Genetic_Diversity-1 ,
                    data = data_timeextinction, family = "poisson")
summary(mod1)
## 
## Call:
## glm(formula = Generation_Eggs ~ Genetic_Diversity - 1, family = "poisson", 
##     data = data_timeextinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.46352  -0.37607   0.00000   0.09066   0.52699  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## Genetic_DiversityMedium   1.5686     0.2041   7.685 1.53e-14 ***
## Genetic_DiversityLow      1.6094     0.1581  10.179  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 103.4310  on 13  degrees of freedom
## Residual deviance:   1.3824  on 11  degrees of freedom
## AIC: 50.301
## 
## Number of Fisher Scoring iterations: 4
datatransformed <- confint(mod1)
exp(datatransformed)
##                            2.5 %   97.5 %
## Genetic_DiversityMedium 3.126628 6.984351
## Genetic_DiversityLow    3.606242 6.713488
exp(1.568)
## [1] 4.797045
exp(1.6094)
## [1] 4.99981

3.3 Survival analysis

str(data_survival)
## 'data.frame':    84 obs. of  6 variables:
##  $ Population       : Factor w/ 84 levels "High1","High13",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Block            : Factor w/ 3 levels "Block3","Block4",..: 1 2 2 2 2 2 2 2 1 2 ...
##  $ Genetic_Diversity: Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Gen_eggs_extinct : int  6 6 6 6 6 6 6 6 6 6 ...
##  $ Survival         : Factor w/ 2 levels "extinct","yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ status           : num  0 0 0 0 0 0 0 0 0 0 ...
#Note: the Surv() function in the {survival} package accepts by default TRUE/FALSE, where TRUE is event and FALSE is censored; 1/0 where 1 is event and 0 is censored; or 2/1 where 2 is event and 1 is censored. Please take care to ensure the event indicator is properly formatted.
  

# The Surv() function from the {survival} package creates a survival object for use as the response in a model formula. There will be one entry for each subject that is the survival time, which is followed by a + if the subject was censored. Let’s look at the first 10 observations:
Surv(data_survival$Gen_eggs_extinct, data_survival$status)
##  [1] 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+
## [26] 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6  6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 5  4  6+ 6+ 5 
## [51] 5  4  6  6+ 6+ 5  6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 5  6+ 4  6+ 6+ 6+ 5  6+ 4 
## [76] 6+ 6+ 6+ 6+ 6  6+ 6+ 6+ 6+
#Surv(data_survival_all$Gen_eggs_extinct, data_survival$status)



# The survfit() function creates survival curves using the Kaplan-Meier method based on a formula. Let’s generate the overall survival curve for the entire cohort, assign it to object s1, and look at the structure using str():
s1 <- survfit(Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival)
str(s1)
## List of 16
##  $ n        : int 84
##  $ time     : num [1:3] 4 5 6
##  $ n.risk   : num [1:3] 84 80 74
##  $ n.event  : num [1:3] 4 6 3
##  $ n.censor : num [1:3] 0 0 71
##  $ surv     : num [1:3] 0.952 0.881 0.845
##  $ std.err  : num [1:3] 0.0244 0.0401 0.0467
##  $ cumhaz   : num [1:3] 0.0476 0.1226 0.1632
##  $ std.chaz : num [1:3] 0.0238 0.0388 0.0453
##  $ type     : chr "right"
##  $ logse    : logi TRUE
##  $ conf.int : num 0.95
##  $ conf.type: chr "log"
##  $ lower    : num [1:3] 0.908 0.814 0.771
##  $ upper    : num [1:3] 0.999 0.953 0.926
##  $ call     : language survfit(formula = Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival)
##  - attr(*, "class")= chr "survfit"
# time: the timepoints at which the curve has a step, i.e. at least one event occurred
# surv: the estimate of survival at the corresponding time


   
# The {ggsurvfit} package works best if you create the survfit object using the included ggsurvfit::survfit2() function, which uses the same syntax to what we saw previously with survival::survfit(). The ggsurvfit::survfit2() tracks the environment from the function call, which allows the plot to have better default values for labeling and p-value reporting.
# 

####
#### PLOT 
####
ggsurvfit::survfit2(Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival) %>% 
  ggsurvfit() +
  ggplot2::labs( x = "Generation",
        y = "Overall survival probability") + 
  add_confidence_interval()  +
  add_risktable()

summary(survfit(Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival))
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     84       4    0.952  0.0232        0.908        0.999
##     5     80       6    0.881  0.0353        0.814        0.953
##     6     74       3    0.845  0.0395        0.771        0.926
##### COMPARISONS ACROSS GROUPS
# We can conduct between-group significance tests using a log-rank test. The log-rank test equally weights observations over the entire follow-up time and is the most common way to compare survival times between groups. There are versions that more heavily weight the early or late follow-up that could be more appropriate depending on the research question (see ?survdiff for different test options).
# 
# We get the log-rank p-value using the survdiff function. For example, we can test whether there was a difference in survival time according to the demography history in the lung data:
  
  

#Comparison 
survdiff(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
## Call:
## survdiff(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, 
##     data = data_survival)
## 
##                           N Observed Expected (O-E)^2/E (O-E)^2/V
## Genetic_Diversity=High   27        0     4.41     4.405      7.00
## Genetic_Diversity=Medium 23        5     3.44     0.707      1.01
## Genetic_Diversity=Low    34        8     5.15     1.571      2.73
## 
##  Chisq= 7  on 2 degrees of freedom, p= 0.03
summary(survfit(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival))
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, 
##     data = data_survival)
## 
##                 Genetic_Diversity=High 
##      time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 
##                 Genetic_Diversity=Medium 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     23       2    0.913  0.0588        0.805        1.000
##     5     21       2    0.826  0.0790        0.685        0.996
##     6     19       1    0.783  0.0860        0.631        0.971
## 
##                 Genetic_Diversity=Low 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     34       2    0.941  0.0404        0.865        1.000
##     5     32       4    0.824  0.0654        0.705        0.962
##     6     28       2    0.765  0.0727        0.635        0.921
####
#### PLOT 
####
plot_survival <- ggsurvfit::survfit2(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival) %>% 
  ggsurvfit(size = 2) +
  labs( x = "Generation",
        y = "Overall survival probability") + 
  add_confidence_interval(alpha = 0.1)  + 
  add_censor_mark(size = 8, shape = "x") +
  scale_fill_manual(name="Demographic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07", "#E7B800")) + 
  scale_color_manual(name="Demographic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07", "#E7B800"))  #+  add_risktable()
plot_survival

# cowplot::save_plot(file = here::here("figures","Extinction_Survival_analysis.pdf"),   
#                    plot_survival, base_height = 10/cm(1),  base_width = 15/cm(1), dpi = 610)


# COX representation 
coxph(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
## Call:
## coxph(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, 
##     data = data_survival)
## 
##                              coef exp(coef)  se(coef)     z     p
## Genetic_DiversityMedium 1.967e+01 3.475e+08 7.239e+03 0.003 0.998
## Genetic_DiversityLow    1.974e+01 3.727e+08 7.239e+03 0.003 0.998
## 
## Likelihood ratio test=11.1  on 2 df, p=0.003883
## n= 84, number of events= 13
coxph(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival) %>% 
  gtsummary::tbl_regression(exp = TRUE) 
Characteristic HR1 95% CI1 p-value
Genetic_Diversity
High — —
Medium 347,496,586 0.00, Inf >0.9
Low 372,735,433 0.00, Inf >0.9
1 HR = Hazard Ratio, CI = Confidence Interval
# The quantity of interest from a Cox regression model is a hazard ratio (HR). The HR represents the ratio of hazards between two groups at any particular point in time. The HR is interpreted as the instantaneous rate of occurrence of the event of interest in those who are still at risk for the event. It is not a risk, though it is commonly mis-interpreted as such. If you have a regression parameter β, then HR = exp(β).
# 
# A HR < 1 indicates reduced hazard of death whereas a HR > 1 indicates an increased hazard of death.
# 
# So the HR = 0.59 implies that 0.59 times as many females are dying as males, at any given time. Stated differently, females have a significantly lower hazard of death than males in these data.







### OTHER REPRESENTATION :  
s2 <- survfit(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
summary(s2)
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, 
##     data = data_survival)
## 
##                 Genetic_Diversity=High 
##      time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 
##                 Genetic_Diversity=Medium 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     23       2    0.913  0.0588        0.805        1.000
##     5     21       2    0.826  0.0790        0.685        0.996
##     6     19       1    0.783  0.0860        0.631        0.971
## 
##                 Genetic_Diversity=Low 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     34       2    0.941  0.0404        0.865        1.000
##     5     32       4    0.824  0.0654        0.705        0.962
##     6     28       2    0.765  0.0727        0.635        0.921
print(s2)
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, 
##     data = data_survival)
## 
##                           n events median 0.95LCL 0.95UCL
## Genetic_Diversity=High   27      0     NA      NA      NA
## Genetic_Diversity=Medium 23      5     NA      NA      NA
## Genetic_Diversity=Low    34      8     NA      NA      NA
# Access to the sort summary table
summary(s2)$table
##                          records n.max n.start events    rmean  se(rmean)
## Genetic_Diversity=High        27    27      27      0 6.000000 0.00000000
## Genetic_Diversity=Medium      23    23      23      5 5.739130 0.12627260
## Genetic_Diversity=Low         34    34      34      8 5.764706 0.09355367
##                          median 0.95LCL 0.95UCL
## Genetic_Diversity=High       NA      NA      NA
## Genetic_Diversity=Medium     NA      NA      NA
## Genetic_Diversity=Low        NA      NA      NA
d <- data.frame(time = s2$time,
                  n.risk = s2$n.risk,
                  n.event = s2$n.event,
                  n.censor = s2$n.censor,
                  surv = s2$surv,
                  upper = s2$upper,
                  lower = s2$lower)
head(d)
##   time n.risk n.event n.censor      surv     upper     lower
## 1    6     27       0       27 1.0000000 1.0000000 1.0000000
## 2    4     23       2        0 0.9130435 1.0000000 0.8048548
## 3    5     21       2        0 0.8260870 0.9964666 0.6848395
## 4    6     19       1       18 0.7826087 0.9707088 0.6309579
## 5    4     34       2        0 0.9411765 1.0000000 0.8653187
## 6    5     32       4        0 0.8235294 0.9621763 0.7048611
plot2 <- survminer::ggsurvplot(s2,
                               pval = TRUE,             # show p-value of log-rank test.
                               #pval.coord  = c(1,0.2), 
                               pval.size = 5, 
                               conf.int = TRUE,         # show confidence intervals for point estimaes of survival curves.
                               # conf.int.style = "step",  # customize style of confidence intervals
                               conf.int.style = "ribbon",  # customize style of confidence intervals 
                               conf.int.alpha = 0.3,  # customize style of confidence intervals
                               censor = TRUE, 
                               censor.shape = "x", 
                               censor.size = 6, 
                               xlab = "Generation",   # customize X axis label.
                               break.time.by = 1,     # break X axis in time intervals by 1
                               ggtheme = theme_light(), # customize plot and risk table with a theme.
                              #risk.table = TRUE, # Add risk table
                              #linetype = "strata", # Change line type by groups
                               risk.table = "abs_pct",  # absolute number and percentage at risk.
                               risk.table.y.text.col = T,# colour risk table text annotations.
                               risk.table.y.text = FALSE,# show bars instead of names in text annotations in legend of risk table.
                               risk.table.col = "strata", # Change risk table color by groups
                               legend.labs = c("Diverse","Intermediate bottleneck","Strong bottleneck"),   
                               legend = c("right"),
                              legend.title = c("Demographic\nhistory:"),   
                               palette = c("#00AFBB", "#E7B800", "#FC4E07")) 



plot2

# 
# cowplot::save_plot(file = here::here("figures","Extinction_Survival_analysis_V2.pdf"),
#                    plot2, base_height = 10/cm(1),  base_width = 15/cm(1), dpi = 610)



# Fit a Cox proportional hazards model
# fit.coxph <- coxph(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
# ggforest(fit.coxph, data = data_survival)

3.4 Population size over time

3.4.1 Polynomial analysis

3.4.1.1 Prelim: test function

####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY WITHOUT GEN 1  
####################### 
###################### REMOVE GEN 1 ##########################33

#If we dont consider Gen=1 
PLOT_Pop_size_average

fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])

fit <- glm(Nb_adults ~ Generation_Eggs, 
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit6 <- glm(Nb_adults~exp(Generation_Eggs), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit7 <- glm(Nb_adults~log(Generation_Eggs), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")





#generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(1,6, length=100)
# plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
#      log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
# lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
# lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
# lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
# lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
# lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
# lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
     data[data$Generation_Eggs>=2,]$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
rsq::rsq(fit,adj=TRUE)
## [1] 0.430959
rsq::rsq(fit2,adj=TRUE)
## [1] 0.5735495
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5723284
rsq::rsq(fit4,adj=TRUE)
## [1] 0.5711286
rsq::rsq(fit5,adj=TRUE)
## [1] 0.5711286
rsq::rsq(fit6,adj=TRUE)
## [1] 0.1285168
rsq::rsq(fit7,adj=TRUE)
## [1] 0.504535
# Best model is fit 2

#data$gen_square <- data$Generation_Eggs^2
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit2 <- glm(Nb_adults~Generation_Eggs + gen_square, 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")

# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
#Add the interaction with Genetic diversity
fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
             gen_square*Genetic_Diversity,
           data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], 
           family = "poisson")


summary(fit2)
## 
## Call:
## glm(formula = Nb_adults ~ Generation_Eggs * Genetic_Diversity + 
##     gen_square * Genetic_Diversity, family = "poisson", data = data[data$Generation_Eggs >= 
##     2 & data$Extinction_finalstatus == "no", ])
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -20.073   -4.181   -0.241    3.315   17.561  
## 
## Coefficients:
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              8.471305   0.053931 157.076  < 2e-16
## Generation_Eggs                         -1.688615   0.031961 -52.833  < 2e-16
## Genetic_DiversityMedium                 -0.593773   0.093819  -6.329 2.47e-10
## Genetic_DiversityLow                     0.157769   0.085778   1.839  0.06588
## gen_square                               0.184959   0.004076  45.383  < 2e-16
## Generation_Eggs:Genetic_DiversityMedium  0.286248   0.056055   5.107 3.28e-07
## Generation_Eggs:Genetic_DiversityLow    -0.241226   0.051392  -4.694 2.68e-06
## Genetic_DiversityMedium:gen_square      -0.050817   0.007248  -7.011 2.36e-12
## Genetic_DiversityLow:gen_square          0.019919   0.006595   3.020  0.00253
##                                            
## (Intercept)                             ***
## Generation_Eggs                         ***
## Genetic_DiversityMedium                 ***
## Genetic_DiversityLow                    .  
## gen_square                              ***
## Generation_Eggs:Genetic_DiversityMedium ***
## Generation_Eggs:Genetic_DiversityLow    ***
## Genetic_DiversityMedium:gen_square      ***
## Genetic_DiversityLow:gen_square         ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 27751  on 352  degrees of freedom
## Residual deviance: 11836  on 344  degrees of freedom
##   (2 observations deleted due to missingness)
## AIC: 14121
## 
## Number of Fisher Scoring iterations: 5
####### Next step: Final step: add genetic diversity 

3.4.1.2 Analysis

#Test oversdispersion
mod0 <- glm(Nb_adults~Generation_Eggs*Genetic_Diversity + 
              gen_square*Genetic_Diversity + 
              Block , 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], 
            family = "poisson")
sqrt(deviance(mod0)/(nobs(mod0)-length(coef(mod0))))
## [1] 5.840179
sigma(mod0)
## [1] 5.840179
# ccl: overdispersion


mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs) , 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], 
            family = "poisson")
blmeco::dispersion_glmer(mod1)
## [1] 1.019726
sigma(mod1)
## [1] 1
#Convergence issue 
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.0111561
mod1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *  
##     Genetic_Diversity + Block + (1 | obs)
##    Data: data[data$Generation_Eggs >= 2 & data$Extinction_finalstatus ==  
##     "no", ]
##       AIC       BIC    logLik  deviance  df.resid 
##  4007.501  4053.899 -1991.751  3983.501       341 
## Random effects:
##  Groups Name        Std.Dev.
##  obs    (Intercept) 0.6643  
## Number of obs: 353, groups:  obs, 353
## Fixed Effects:
##                             (Intercept)  
##                                 9.22506  
##                         Generation_Eggs  
##                                -2.10047  
##                 Genetic_DiversityMedium  
##                                -0.65107  
##                    Genetic_DiversityLow  
##                                 0.10736  
##                              gen_square  
##                                 0.23715  
##                             BlockBlock4  
##                                -0.12145  
##                             BlockBlock5  
##                                -0.48597  
## Generation_Eggs:Genetic_DiversityMedium  
##                                 0.32468  
##    Generation_Eggs:Genetic_DiversityLow  
##                                -0.24262  
##      Genetic_DiversityMedium:gen_square  
##                                -0.06201  
##         Genetic_DiversityLow:gen_square  
##                                 0.01594  
## convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#Test add Population random
mod0 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs) + (1|Population), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], 
            family = "poisson")

max(abs(with(mod0@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.006369406
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.0111561
anova(mod0, mod1, test ="Chisq")  
## Data: data[data$Generation_Eggs >= 2 & data$Extinction_finalstatus ==  ...
## Models:
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## mod1:     Genetic_Diversity + Block + (1 | obs)
## mod0: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## mod0:     Genetic_Diversity + Block + (1 | obs) + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1   12 4007.5 4053.9 -1991.8   3983.5                         
## mod0   13 3993.9 4044.2 -1984.0   3967.9 15.611  1  7.782e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
########## LRT and Final model: 
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs)  + (1|Population), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + Genetic_Diversity + 
                Block  + (1|obs)  + (1|Population), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod1, mod2, test ="Chisq")  
## Data: data[data$Generation_Eggs >= 2 & data$Extinction_finalstatus ==  ...
## Models:
## mod2: Nb_adults ~ Generation_Eggs + gen_square + Genetic_Diversity + 
## mod2:     Block + (1 | obs) + (1 | Population)
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## mod1:     Genetic_Diversity + Block + (1 | obs) + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod2    9 3997.6 4032.4 -1989.8   3979.6                       
## mod1   13 3993.9 4044.2 -1984.0   3967.9 11.663  4    0.02004 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(2,6, length = 100),each=3),
                       Block = "Block4",
                       gen_square = (rep(seq(2,6, length = 100),each=3))^2,
                       Estimates = NA)

filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

emmeans::emmeans(mod1, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean     SE  df asymp.LCL asymp.UCL
##  High                4.89 0.0763 Inf      4.71      5.08
##  Medium              4.42 0.0930 Inf      4.20      4.65
##  Low                 4.31 0.0801 Inf      4.12      4.51
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.471 0.120 Inf 3.917   0.0003 
##  High - Low       0.581 0.110 Inf 5.287   <.0001 
##  Medium - Low     0.109 0.122 Inf 0.899   0.6408 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.4.1.3 Extra: with gen1

###################################################################
###################################################################
###################################################################
###################################################################
####################  EXTRA: INCLUDING GEN 1 ###########################
PLOT_Pop_size_average

fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")


#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
     data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
summary(fit)$adj.r.squared
## NULL
summary(fit2)$adj.r.squared
## NULL
summary(fit3)$adj.r.squared
## NULL
summary(fit4)$adj.r.squared # best r square
## NULL
summary(fit5)$adj.r.squared 
## NULL
summary(fit6)$adj.r.squared 
## NULL
summary(fit7)$adj.r.squared 
## NULL
rsq::rsq(fit,adj=TRUE)
## [1] 0.1076945
rsq::rsq(fit2,adj=TRUE)
## [1] 0.1498895
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5256167
rsq::rsq(fit4,adj=TRUE)
## [1] 0.5982229
rsq::rsq(fit5,adj=TRUE)
## [1] 0.5992474
rsq::rsq(fit6,adj=TRUE)
## [1] 0.07238302
rsq::rsq(fit7,adj=TRUE)
## [1] 0.04993927
#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")


####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY and GEN 1 
####################### 

#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)
## [1] 5.776877
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)
## [1] 1.086207
#Convergence issue 
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.05822911
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(1,6, length = 100),each=3),
                       Block = "Block4", 
                       Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity + 
#              gen_square*Genetic_Diversity, 
#            data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity + 
#              gen_square, 
#            data = data[data$Generation_Eggs>=2,])

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    9 5628.3 5666.1 -2805.2   5610.3                         
## mod    17 5611.8 5683.0 -2788.9   5577.8 32.578  8  7.336e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Signifcant interaction 
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5611.8   5683.0  -2788.9   5577.8      471 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.01761 -0.03977  0.00658  0.05613  0.26106 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.66352  0.8146  
##  Block  (Intercept) 0.07891  0.2809  
## Number of obs: 488, groups:  obs, 488; Block, 3
## 
## Fixed effects:
##                                                                       Estimate
## (Intercept)                                                          -1.993366
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.860627
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.148606
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.942526
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.059012
## Genetic_DiversityMedium                                               0.024770
## Genetic_DiversityLow                                                 -0.788813
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  0.061583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.079133
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.003441
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.001629
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     1.488961
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    -0.832919
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.138298
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    -0.006970
##                                                                      Std. Error
## (Intercept)                                                            1.252012
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           2.008686
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                           1.044443
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           0.215778
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                           0.015351
## Genetic_DiversityMedium                                                1.849227
## Genetic_DiversityLow                                                   1.649489
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   3.000571
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   1.565125
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.324403
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.023149
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      2.669129
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      1.389632
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.287697
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.020517
##                                                                      z value
## (Intercept)                                                           -1.592
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           5.407
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -4.930
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           4.368
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -3.844
## Genetic_DiversityMedium                                                0.013
## Genetic_DiversityLow                                                  -0.478
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.021
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  -0.051
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  -0.011
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.070
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.558
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     -0.599
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.481
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     -0.340
##                                                                      Pr(>|z|)
## (Intercept)                                                          0.111355
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         6.41e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         8.24e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         1.25e-05
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         0.000121
## Genetic_DiversityMedium                                              0.989313
## Genetic_DiversityLow                                                 0.632496
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.983625
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.959676
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.991538
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.943917
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    0.576950
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    0.548919
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    0.630726
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    0.734088
##                                                                         
## (Intercept)                                                             
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 4.38011 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
summary(mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5628.3   5666.1  -2805.2   5610.3      479 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.03098 -0.04847  0.01332  0.06982  0.21353 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.71155  0.8435  
##  Block  (Intercept) 0.07597  0.2756  
## Number of obs: 488, groups:  obs, 488; Block, 3
## 
## Fixed effects:
##                                              Estimate Std. Error z value
## (Intercept)                                  -1.90446    0.82550  -2.307
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.54204    1.32511   8.710
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.54378    0.69569  -7.969
## stats::poly(Generation_Eggs, 4, raw = TRUE)3  1.00561    0.14464   6.952
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.06197    0.01033  -5.999
## Genetic_DiversityMedium                      -0.58397    0.10123  -5.769
## Genetic_DiversityLow                         -0.67700    0.09165  -7.387
##                                              Pr(>|z|)    
## (Intercept)                                    0.0211 *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 1.60e-15 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 3.59e-12 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 1.99e-09 ***
## Genetic_DiversityMedium                      7.98e-09 ***
## Genetic_DiversityLow                         1.50e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                    (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.966                                      
## s::(G_E,4,r=TRUE)2  0.942 -0.992                               
## s::(G_E,4,r=TRUE)3 -0.915  0.976             -0.995            
## s::(G_E,4,r=TRUE)4  0.889 -0.956              0.984            
## Gntc_DvrstM        -0.060  0.004             -0.005            
## Gntc_DvrstL        -0.065  0.004             -0.005            
##                    s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1                                             
## s::(G_E,4,r=TRUE)2                                             
## s::(G_E,4,r=TRUE)3                                             
## s::(G_E,4,r=TRUE)4 -0.997                                      
## Gntc_DvrstM         0.005             -0.005                   
## Gntc_DvrstL         0.005             -0.006              0.495
## convergence code: 0
## Model failed to converge with max|grad| = 0.0615683 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.55 0.206 Inf      4.06      5.04
##  Medium              3.94 0.217 Inf      3.42      4.46
##  Low                 3.68 0.201 Inf      3.20      4.16
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.609 0.189 Inf 3.221   0.0037 
##  High - Low       0.866 0.170 Inf 5.078   <.0001 
##  Medium - Low     0.257 0.186 Inf 1.379   0.3517 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
####################### 
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
####################### 
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
            data = data[data$Extinction_finalstatus=="no",], family = "poisson")

filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")

PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates_Withoutextpop, 
                               colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
  geom_point(data = data[data$Extinction_finalstatus=="no",], 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), 
             data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), 
              data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod1    9 4759.9 4796.3 -2370.9   4741.9                       
## mod    17 4756.3 4825.1 -2361.1   4722.3 19.599  8    0.01196 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data[data$Extinction_finalstatus == "no", ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   4756.3   4825.1  -2361.1   4722.3      407 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57063 -0.04879  0.00393  0.06967  0.29650 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.36624  0.6052  
##  Block  (Intercept) 0.02456  0.1567  
## Number of obs: 424, groups:  obs, 424; Block, 3
## 
## Fixed effects:
##                                                                        Estimate
## (Intercept)                                                          -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.0583188
## Genetic_DiversityMedium                                               1.1892488
## Genetic_DiversityLow                                                  0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0006044
##                                                                      Std. Error
## (Intercept)                                                           0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                          1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          0.0120893
## Genetic_DiversityMedium                                               1.5049889
## Genetic_DiversityLow                                                  1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0173582
##                                                                      z value
## (Intercept)                                                           -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -4.824
## Genetic_DiversityMedium                                                0.790
## Genetic_DiversityLow                                                   0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.035
##                                                                      Pr(>|z|)
## (Intercept)                                                            0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         1.41e-06
## Genetic_DiversityMedium                                                0.4294
## Genetic_DiversityLow                                                   0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.9722
##                                                                         
## (Intercept)                                                          *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.53 0.133 Inf      4.21      4.85
##  Medium              4.30 0.151 Inf      3.94      4.66
##  Low                 4.01 0.137 Inf      3.68      4.33
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.230 0.153 Inf 1.507   0.2874 
##  High - Low       0.523 0.140 Inf 3.746   0.0005 
##  Medium - Low     0.293 0.157 Inf 1.861   0.1503 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.4.1.4 Extra: Pop. size G6

############ 
###### Analysis
############ 
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),  
#       data = data[data$Generation_Eggs==6,], family = "poisson")


mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])

mod1 <- lm(log(Nb_adults) ~ Block , 
             data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])

anova(mod0, mod1) # Effect of genetic diversity 
## Analysis of Variance Table
## 
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     66 29.706                                  
## 2     68 42.118 -2   -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])


summary(mod0)
## 
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data[data$Generation_Eggs == 
##     6 & data$Extinction_finalstatus == "no", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2417 -0.2947  0.1088  0.4381  1.2624 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              5.11532    0.17980  28.451  < 2e-16 ***
## Genetic_DiversityMedium -0.94813    0.20540  -4.616 1.86e-05 ***
## Genetic_DiversityLow    -0.80532    0.18719  -4.302 5.72e-05 ***
## BlockBlock4             -0.02077    0.18447  -0.113   0.9107    
## BlockBlock5             -0.53922    0.21414  -2.518   0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared:  0.3362, Adjusted R-squared:  0.2959 
## F-statistic: 8.356 on 4 and 66 DF,  p-value: 1.623e-05
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE df lower.CL upper.CL
##  High                4.93 0.130 66     4.61     5.25
##  Medium              3.98 0.159 66     3.59     4.37
##  Low                 4.12 0.136 66     3.79     4.46
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE df t.ratio p.value
##  High - Medium    0.948 0.205 66  4.616  0.0001 
##  High - Low       0.805 0.187 66  4.302  0.0002 
##  Medium - Low    -0.143 0.207 66 -0.690  0.7704 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.4.2 GAMMM Model

3.4.2.1 Basic model

# https://jacolienvanrij.com/Tutorials/GAMM.html
# install.packages("itsadug")
# packageVersion('mgcv')
# packageVersion('itsadug')

data_dyn <- data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]
str(data_dyn)
## 'data.frame':    355 obs. of  33 variables:
##  $ Population            : Factor w/ 84 levels "High1","High13",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ Block                 : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ Old_Label             : chr  "B3 All red mix " "B3 All red mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "7/7 BE B3 | G5 High 1" "6/2 BE B3 | G4 High 1" "4/28 BE B3 | G3 High 1" "BE B3 3/24 | G2 High1" ...
##  $ Genetic_Diversity     : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Generation_Parents    : int  4 3 2 1 5 5 3 2 1 4 ...
##  $ Generation_Eggs       : int  5 4 3 2 6 6 4 3 2 5 ...
##  $ Date_Census           : chr  "7/7/21" "6/2/21" "4/28/21" "3/24/21" ...
##  $ Date_Start_Eggs       : chr  "7/7/21" "6/2/21" "4/28/21" "3/24/21" ...
##  $ Date_End_Eggs         : chr  "7/8/21" "6/3/21" "4/29/21" "3/25/21" ...
##  $ Image_Box1            : chr  "DSC_0779" "DSC_0206" "DSC_0585" "DSC_0920" ...
##  $ Image_Box2            : chr  "DSC_0780" "DSC_0207" "DSC_0586" "DSC_0921" ...
##  $ MC_Box1               : num  26.6 85.4 82.3 266.9 100.1 ...
##  $ MC_Box2               : num  5.15 90.72 78.01 149.8 74.38 ...
##  $ MC_Total_Adults       : num  31.8 176.1 160.3 416.7 174.5 ...
##  $ HC_Box1               : int  27 79 83 253 64 47 71 60 101 47 ...
##  $ HC_Box2               : int  4 96 82 135 58 80 118 113 214 44 ...
##  $ HC_Total_Adults       : int  31 175 165 388 122 127 189 173 315 91 ...
##  $ Nb_parents            : num  175 165 388 100 31 91 173 315 100 189 ...
##  $ Growth_rate           : num  0.177 1.061 0.425 3.88 3.935 ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  31 175 165 388 122 127 189 173 315 91 ...
##  $ Lambda                : num  0.177 1.061 0.425 3.88 3.935 ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 337 253 169 85 421 447 279 195 111 363 ...
##  $ gen_square            : num  25 16 9 4 36 36 16 9 4 25 ...
##  $ He_start              : num  0.999 0.999 0.999 0.999 0.999 ...
##  $ He_lost_gen_t         : num  0.984 0.997 0.997 0.999 0.996 ...
##  $ He_remain             : num  0.971 0.986 0.989 0.992 0.967 ...
##  $ He_exp                : num  0.968 0.968 0.968 0.968 0.968 ...
##  $ He_end                : num  0.967 0.967 0.967 0.967 0.967 ...
##  $ ID_extinction         : chr  "extant" "extant" "extant" "extant" ...
# 1-dimensional smooth of Time:
#m1 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs), data=data_dyn)
#Meaning the GAM has failed to construct a sensible smooth term for this data. You can limit the maximum df of the smooth using the k parameter:
m1 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, k=3), data=data_dyn)

#number of ‘knots’. This parameter determines the upper bound of the number of underlying base functions being used to build up the curve. Thus, this parameter constraints the wigglyness of a smooth, or - as a metaphor - the number of bowpoints of a curve.
#Note that the model will base the number of base functions (reflected in the edf of the summary) on the data with the setting for k as upper bound. By default, the value of k for s() is around 9, and for te() and ti() 5 per dimension. Importantly, the value of k should be at most one less than the number of unique data points, otherwise it will fit the density of that predictor.

# interaction with non-isotropicterms:
#m2 <- gam(Nb_adults ~ te(Generation_Eggs, Genetic_Diversity), data=data_dyn)

# decompose the same interaction in main effects + interaction:
# (note: include main effects)
# m3 <- gam(Nb_adults ~ s(Generation_Eggs) + s(Genetic_Diversity) +
#           + ti(Generation_Eggs, Genetic_Diversity), 
#           data=data_dyn)

# s() : for modeling a 1-dimensional smooth, or for modeling isotropic interactions (variables are measured in same units and on same scale)
# te(): for modeling 2- or n-dimensional interaction surfaces of variables that are not isotropic (but see info about d parameter below). Includes ‘main’ effects.
# ti(): for modeling 2- or n-dimensional interaction surfaces that do not include the ‘main effects

########### 
########### ADD FACTOR EFFECT: 
########### 

# How to set up the interaction depends on the type of grouping predictor:
# with factor include intercept difference: Group + s(Time, by=Group).
# with ordered factor include intercept difference and reference smooth: Group + s(Time) + s(Time, by=Group).
# with binary predictor include reference smooth: s(Time) + s(Time, by=IsGroupChildren).

m.factor <- mgcv::gam(Nb_adults ~ Genetic_Diversity +  s(Generation_Eggs, k=3) + 
            s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
summary(m.factor)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Nb_adults ~ Genetic_Diversity + s(Generation_Eggs, k = 3) + s(Generation_Eggs, 
##     by = Genetic_Diversity, k = 3)
## 
## Parametric coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              172.606      5.715  30.205  < 2e-16 ***
## Genetic_DiversityMedium  -46.636      9.015  -5.173 3.92e-07 ***
## Genetic_DiversityLow     -55.236      8.160  -6.770 5.56e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                              edf Ref.df       F  p-value    
## s(Generation_Eggs)                         1.714  1.746 129.456  < 2e-16 ***
## s(Generation_Eggs):Genetic_DiversityHigh   1.467  1.667   4.035 0.013248 *  
## s(Generation_Eggs):Genetic_DiversityMedium 0.750  0.750  14.965 0.000895 ***
## s(Generation_Eggs):Genetic_DiversityLow    0.750  0.750   8.142 0.013939 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 10/11
## R-sq.(adj) =  0.615   Deviance explained = 62.2%
## GCV = 4473.2  Scale est. = 4375.8    n = 353
#The smooth terms summary shows three smooths that are significantly different from 0. 
#We cannot conclude from the summary that the two curves are different from each other.


m.factor <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=3) + 
            s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
summary(m.factor)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 3) + 
##     s(Generation_Eggs, by = Genetic_Diversity, k = 3)
## 
## Parametric coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              181.701      7.897  23.008  < 2e-16 ***
## Genetic_DiversityMedium  -47.619      9.023  -5.277 2.33e-07 ***
## Genetic_DiversityLow     -58.396      8.238  -7.088 7.77e-12 ***
## BlockBlock4               -6.052      8.108  -0.746   0.4559    
## BlockBlock5              -21.941      9.421  -2.329   0.0204 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                              edf Ref.df       F  p-value    
## s(Generation_Eggs)                         1.714  1.746 130.723  < 2e-16 ***
## s(Generation_Eggs):Genetic_DiversityHigh   1.472  1.670   4.118 0.012347 *  
## s(Generation_Eggs):Genetic_DiversityMedium 0.750  0.750  15.121 0.000842 ***
## s(Generation_Eggs):Genetic_DiversityLow    0.750  0.750   8.228 0.013454 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 12/13
## R-sq.(adj) =  0.619   Deviance explained = 62.8%
## GCV = 4452.8  Scale est. = 4330.7    n = 353
########### 
########### ADD RANDOM EFFECT: 
########### 

# random intercepts adjust the height of other modelterms with a constant value: s(Subject, bs="re").
# random slopes adjust the slope ofthe trend of a numeric predictor: s(Subject, Time, bs="re").
# random smooths adjust the trend of a numeric predictor in a nonlinear way: s(Time, Subject, bs="fs", m=1).

m.mixed <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) + 
            s(Generation_Eggs, by= Genetic_Diversity, k=3) + 
             s(Population, bs="re") , data=data_dyn)
summary(m.mixed)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 3) + 
##     s(Generation_Eggs, by = Genetic_Diversity, k = 3) + s(Population, 
##     bs = "re")
## 
## Parametric coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              181.781      8.979  20.245  < 2e-16 ***
## Genetic_DiversityMedium  -47.725     10.259  -4.652 4.79e-06 ***
## Genetic_DiversityLow     -58.512      9.362  -6.250 1.29e-09 ***
## BlockBlock4               -6.094      9.217  -0.661   0.5089    
## BlockBlock5              -21.800     10.707  -2.036   0.0426 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                               edf Ref.df       F  p-value    
## s(Generation_Eggs)                          1.714  1.746 140.096  < 2e-16 ***
## s(Generation_Eggs):Genetic_DiversityHigh    1.492  1.681   4.504 0.008956 ** 
## s(Generation_Eggs):Genetic_DiversityMedium  0.750  0.750  16.246 0.000547 ***
## s(Generation_Eggs):Genetic_DiversityLow     0.750  0.750   8.868 0.010341 *  
## s(Population)                              18.456 66.000   0.389 0.032725 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 83/84
## R-sq.(adj) =  0.646   Deviance explained = 67.3%
## GCV = 4377.5  Scale est. = 4028.3    n = 353

3.4.2.2 Model selection

########### 
########### AIC
########### 
m1 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
m2 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=4), data=data_dyn)
m3 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=5), data=data_dyn)
m4 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
m5 <- mgcv::gam(Nb_adults ~ Block+ 
                  s(Generation_Eggs, by= Genetic_Diversity, k=4), data=data_dyn)
m6 <- mgcv::gam(Nb_adults ~ Block +  s(Generation_Eggs, k=3) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5), data=data_dyn)
m7 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=3) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5), data=data_dyn)
m8 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=3) + 
                   s(Population, bs="re"), data=data_dyn)
m9 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=4) + 
                   s(Population, bs="re"), data=data_dyn)
m10 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=5)+
                  s(Population, bs="re"),data=data_dyn)
m11 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=3)+
                    s(Population, bs="re"), data=data_dyn)
m12 <- mgcv::gam(Nb_adults ~ Block +  
                  s(Generation_Eggs, by= Genetic_Diversity, k=4) + 
                    s(Population, bs="re"), data=data_dyn)
m13 <- mgcv::gam(Nb_adults ~ Block +  s(Generation_Eggs, k=3) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
m14 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=3) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
m15 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=5) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
m16 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=4) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=4) + 
                    s(Population, bs="re"), data=data_dyn)
m17 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=3) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=3) + 
                    s(Population, bs="re"), data=data_dyn)

AIC(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17)
##            df      AIC
## m1   7.936534 4017.558
## m2  10.081634 4011.546
## m3  11.048573 4011.346
## m4  11.944931 3970.451
## m5  12.082285 4013.676
## m6  11.674632 4012.361
## m7  13.980856 3961.220
## m8  46.262085 3975.939
## m9  49.630906 3965.567
## m10 51.228589 3963.659
## m11 30.598171 3962.136
## m12 50.983869 3966.744
## m13 51.608090 3964.162
## m14 35.378282 3950.158
## m15 33.635622 3947.089
## m16 32.497862 3948.671
## m17 29.161907 3961.031
#https://www.seascapemodels.org/rstats/2021/03/27/common-GAM-problems.html




# Real comparison 
mod0 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=5) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
mod1 <- mgcv::gam(Nb_adults ~ Block +  s(Generation_Eggs, k=5) + 
                  s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
mod2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
mod3 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  
                    s(Generation_Eggs, by= Genetic_Diversity, k=3) + 
                    s(Population, bs="re"), data=data_dyn)
mod4 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  
                    s(Generation_Eggs, by= Genetic_Diversity, k=4) + 
                    s(Population, bs="re"), data=data_dyn)
mod5 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  
                    s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
mod7 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +   s(Generation_Eggs, k=3) +
                    s(Generation_Eggs, by= Genetic_Diversity, k=3) + 
                    s(Population, bs="re"), data=data_dyn)
mod8 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +   s(Generation_Eggs, k=4) +
                    s(Generation_Eggs, by= Genetic_Diversity, k=4) + 
                    s(Population, bs="re"), data=data_dyn)
mod9 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +   s(Generation_Eggs, k=5) +
                    s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)

AIC(mod1, mod2, mod3, mod4, mod5, mod7, mod8, mod9)
##            df      AIC
## mod1 49.81430 3961.188
## mod2 30.46099 3947.366
## mod3 30.59817 3962.136
## mod4 34.56955 3952.491
## mod5 36.35822 3950.893
## mod7 29.16191 3961.031
## mod8 32.49786 3948.671
## mod9 33.63562 3947.089
mod2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +  s(Generation_Eggs, k=5) + 
                    s(Population, bs="re"), data=data_dyn)
mod9 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +   s(Generation_Eggs, k=5) +
                    s(Generation_Eggs, by= Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), data=data_dyn)


AIC(mod2, mod8, mod9)
##            df      AIC
## mod2 30.46099 3947.366
## mod8 32.49786 3948.671
## mod9 33.63562 3947.089
visreg::visreg(mod9, xvar = "Generation_Eggs",
       by = "Genetic_Diversity", data = data_dyn,
       method = "REML")

visreg::visreg(mod2, xvar = "Generation_Eggs",
       by = "Genetic_Diversity", data = data_dyn,
       method = "REML")

3.4.2.3 Testing for significance

# Three important methods to determine the significance of predictors:
## -1 Model-comparison procedure (using function compareML).
## -2 Test statistics of the smooth terms as presented in the summary.
## -3 Visualization of the smooth terms (e.g., difference plots and difference surfaces) but useful only with significative interaction https://jacolienvanrij.com/Tutorials/GAMM.html 


###### 1- Model comparison 

mod_1 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +   s(Generation_Eggs, k=5) +
                    s(Generation_Eggs, by = Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), method="ML", data=data_dyn)
mod_2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block  +  s(Generation_Eggs, k=5) +
                    s(Population, bs="re"), data=data_dyn, method='ML')

itsadug::compareML(mod_1, mod_2)
## mod_1: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Generation_Eggs, by = Genetic_Diversity, k = 5) + s(Population, 
##     bs = "re")
## 
## mod_2: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Population, bs = "re")
## 
## Chi-square test of ML scores
## -----
##   Model    Score Edf Difference    Df p.value Sig.
## 1 mod_2 1975.142   8                              
## 2 mod_1 1974.369  14      0.772 6.000   0.957     
## 
## AIC difference: 2.51, model mod_2 has lower AIC.
AIC(mod_1, mod_2)
##             df      AIC
## mod_1 34.90059 3956.820
## mod_2 30.31497 3954.313
mod_2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block  +  s(Generation_Eggs, k=5) +
                    s(Population, bs="re"), data=data_dyn, method='ML')
mod_3 <- mgcv::gam(Nb_adults ~ Block  +  s(Generation_Eggs, k=5) +
                    s(Population, bs="re"), data=data_dyn, method='ML')
itsadug::compareML(mod_3, mod_2)
## mod_3: Nb_adults ~ Block + s(Generation_Eggs, k = 5) + s(Population, 
##     bs = "re")
## 
## mod_2: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Population, bs = "re")
## 
## Chi-square test of ML scores
## -----
##   Model    Score Edf Difference    Df   p.value Sig.
## 1 mod_3 1993.015   6                                
## 2 mod_2 1975.142   8     17.873 2.000 1.729e-08  ***
## 
## AIC difference: 10.73, model mod_2 has lower AIC.
mod_final <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block  +  
                         s(Generation_Eggs, k=5) +
                         s(Population, bs="re"), data=data_dyn, method='ML')


mod_final <- mgcv::gam(Nb_adults ~ Genetic_Diversity-1 + Block  +  
                         s(Generation_Eggs, k=5) +
                         s(Population, bs="re"), data=data_dyn, method='ML')


###### 2- Summary of the model
# With factors such as Group we cannot use the summary to test differences between groups. However, to investigate differences between groups, we could use difference smooths with ordered factors or binary predictors.
data_dyn$OFGroup_diversity <- as.factor(data_dyn$Genetic_Diversity)
# change factor to ordered factor:
data_dyn$OFGroup_diversity <- as.ordered(data_dyn$OFGroup)
# change contrast to treatment coding (difference curves)
contrasts(data_dyn$OFGroup_diversity) <- 'contr.treatment'
# Inspect contrasts:
contrasts(data_dyn$OFGroup_diversity)
##        Medium Low
## High        0   0
## Medium      1   0
## Low         0   1
mod_1 <- mgcv::gam(Nb_adults ~ OFGroup_diversity + Block +   s(Generation_Eggs, k=5) +
                    s(Generation_Eggs, by = OFGroup_diversity, k=5) + 
                    s(Population, bs="re"), method="ML", data=data_dyn)


summary(mod_1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Nb_adults ~ OFGroup_diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Generation_Eggs, by = OFGroup_diversity, k = 5) + s(Population, 
##     bs = "re")
## 
## Parametric coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              181.737      8.773  20.716  < 2e-16 ***
## OFGroup_diversityMedium  -47.779     10.023  -4.767 2.83e-06 ***
## OFGroup_diversityLow     -58.470      9.148  -6.392 5.70e-10 ***
## BlockBlock4               -6.045      9.005  -0.671   0.5025    
## BlockBlock5              -21.621     10.462  -2.067   0.0396 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                               edf Ref.df      F p-value    
## s(Generation_Eggs)                          3.755  3.963 92.191  <2e-16 ***
## s(Generation_Eggs):OFGroup_diversityMedium  1.217  1.398  0.130  0.6960    
## s(Generation_Eggs):OFGroup_diversityLow     1.001  1.001  0.226  0.6353    
## s(Population)                              17.712 66.000  0.387  0.0208 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.656   Deviance explained = 68.3%
## -ML = 1974.8  Scale est. = 3905.5    n = 353
#We confirm: No interaction significative
# The summary now provides the significance of the difference terms: Medium and High does not show a significantly different trend over Time (F=0.130, edf=1.21, p=0.69), Low and High does not showa significantly different trend over Time (F=0.226, edf=1.001, p=0.02). 

mod_2 <- mgcv::gam(Nb_adults ~ OFGroup_diversity + Block +   s(Generation_Eggs, k=5) +
                    s(Population, bs="re"), method="ML", data=data_dyn)
summary(mod_2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Nb_adults ~ OFGroup_diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Population, bs = "re")
## 
## Parametric coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              181.743      8.751  20.768  < 2e-16 ***
## OFGroup_diversityMedium  -47.779      9.999  -4.778 2.67e-06 ***
## OFGroup_diversityLow     -58.483      9.125  -6.409 5.12e-10 ***
## BlockBlock4               -6.060      8.983  -0.675   0.5004    
## BlockBlock5              -21.617     10.436  -2.071   0.0391 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df       F p-value    
## s(Generation_Eggs)  3.752  3.963 147.751  <2e-16 ***
## s(Population)      17.500 66.000   0.383  0.0205 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.656   Deviance explained = 68.1%
## -ML = 1975.1  Scale est. = 3903.6    n = 353
###### 2-  Visualization of the smooth terms
# seful only with significative interaction 
# Which is not the case here
# But see: https://jacolienvanrij.com/Tutorials/GAMM.html

3.4.2.4 Plot

data_dyn <- data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]

mod_final <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + 
                     s(Generation_Eggs, k=5) +  
                     s(Population, bs="re"), data=data_dyn)


# Plot  block effect
mod_block_P <- tidymv::predict_gam(mod_final, 
                               exclude_terms = list("s(Population)"), 
                               values = list(Population = NULL))




ggplot(data = mod_block_P, aes(Generation_Eggs, fit)) +
    facet_wrap(~ Block) +
    tidymv::geom_smooth_ci(Genetic_Diversity)

# Plot Diversity effect
mod_A_p <- tidymv::predict_gam(mod_final, 
                               exclude_terms = list("s(Population)", "Block"), 
                               values = list(Population = NULL, Block = NULL))

ggplot(data = mod_A_p, aes(Generation_Eggs, fit)) +
   tidymv::geom_smooth_ci(Genetic_Diversity) 

## Prediction for main text 
#population size at the end
as.data.frame(mod_A_p[mod_A_p$Generation_Eggs==6,])
##   Genetic_Diversity  Block Generation_Eggs Population       fit   se.fit
## 1              High Block4               6      High1 139.87617 10.27320
## 2            Medium Block4               6      High1  92.07954 11.62504
## 3               Low Block4               6      High1  81.37454 10.46044
#minimum population size 
tempdata <- as.data.frame(mod_A_p)
tapply(mod_A_p$fit,mod_A_p$Genetic_Diversity,min)
##      High    Medium       Low 
## 113.09324  65.29661  54.59161
tempdata[tempdata$Genetic_Diversity=="High"&tempdata$fit<=114,]
##     Genetic_Diversity  Block Generation_Eggs Population      fit   se.fit
## 91               High Block4        4.448980      High1 113.5282 9.489786
## 94               High Block4        4.530612      High1 113.2052 9.470716
## 97               High Block4        4.612245      High1 113.0932 9.515032
## 100              High Block4        4.693878      High1 113.1954 9.615687
## 103              High Block4        4.775510      High1 113.5150 9.753966
## PLOT

plot_GAM_Ucurve <- ggplot(data = mod_A_p, aes(Generation_Eggs, fit)) +
  geom_point(data = data[data$Extinction_finalstatus=="no",], 
             aes(x = Generation_Eggs,
                                y = Nb_adults,
                               group = Genetic_Diversity,
                               # group = Population,
                                colour = Genetic_Diversity),
             position = position_dodge2(0.3), 
             size = 0.6, alpha = 0.9) +
  tidymv::geom_smooth_ci(Genetic_Diversity, 
                 size = 1.5, linetype = 1)  + 
  ylab("Population size") +
  xlab("Generation") +
  scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"),
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
                   # values = c("#685338", "#FEA698", "#7BC7AD")) +
                   # values = c( "#51392C","#DA5837", "#4C9E97")) +
                   # values = c( "#00AFBB", "#E7B800", "#FC4E07")) +
  
                   #  values = c("#9AD3CA", "#929FD1", "#E7C7D7")) +
                   #values = c("#227C9D", "#FFCB77", "#FE6D73")) +
                     values = c("#00AFBB","#FC4E07","#E7B800")) +
 theme_LO 


plot_GAM_Ucurve

# 
# cowplot::save_plot(file = here::here("figures","PopulationSize_GAM.pdf"),
#                    plot_GAM_Ucurve, base_height = 10/cm(1),           
#                    base_width = 15/cm(1), dpi = 610)
# 

  

#To estimate confidence interval at the end of the experiment
#Default on tidymv::geom_smooth_ci is ci_z = 1.96
#population size at the end
data_gam_end <- as.data.frame(mod_A_p[mod_A_p$Generation_Eggs==6,])
data_gam_end$CIlower <- data_gam_end$fit-(1.96*data_gam_end$se.fit)
data_gam_end$CIupper <- data_gam_end$fit+(1.96*data_gam_end$se.fit)
data_gam_end
##   Genetic_Diversity  Block Generation_Eggs Population       fit   se.fit
## 1              High Block4               6      High1 139.87617 10.27320
## 2            Medium Block4               6      High1  92.07954 11.62504
## 3               Low Block4               6      High1  81.37454 10.46044
##     CIlower  CIupper
## 1 119.74070 160.0116
## 2  69.29445 114.8646
## 3  60.87207 101.8770
#Confirmation on the plot: 
# plot_GAM_Ucurve + 
#   ylim (50,170) + 
#   xlim (5.9,6.1)

3.5 Variance over time

3.5.1 Residuals residuals

#data$gen_square <- data$Generation_Eggs^2
data_G2_all <- data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]
data_G2_all <- data_G2_all[complete.cases(data_G2_all$Nb_adults), ]

Initial_model<- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data_G2_all, 
            family = "poisson")
data_G2_all$resid <- resid(Initial_model)

PLOT_Pop_size_predict

#PLOT RESIDUALS

ggplot2::ggplot(data_G2_all, aes(x = factor(Generation_Eggs),  
                                 y = resid, 
                                 colour = Genetic_Diversity)) + 
  geom_boxplot(outlier.color = NA) +
  geom_point(position = position_dodge(0.75), size = 0.7, alpha = 0.5) +
  scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#E7B800","#FC4E07")) +
  ylab("Resid(pop size)") +
  xlab("Generation") +
  theme_LO

########### 
########### OPTION A: Test random effects 
########### 

# Test variances with residuals and random effects: 
model_withrandomeffect <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs) + (1|Genetic_Diversity:Generation_Eggs) , 
            data = data_G2_all, 
            family = "poisson")
model_withoutrandomeffect <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs) , 
            data = data_G2_all, 
            family = "poisson")

anova(model_withrandomeffect,model_withoutrandomeffect, test = "Chisq")
## Data: data_G2_all
## Models:
## model_withoutrandomeffect: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## model_withoutrandomeffect:     Genetic_Diversity + Block + (1 | obs)
## model_withrandomeffect: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## model_withrandomeffect:     Genetic_Diversity + Block + (1 | obs) + (1 | Genetic_Diversity:Generation_Eggs)
##                           npar    AIC    BIC  logLik deviance Chisq Df
## model_withoutrandomeffect   12 4007.5 4053.9 -1991.8   3983.5         
## model_withrandomeffect      13 4009.5 4059.8 -1991.8   3983.5 3e-04  1
##                           Pr(>Chisq)
## model_withoutrandomeffect           
## model_withrandomeffect        0.9859
AIC(model_withrandomeffect)
## [1] 4009.501
AIC(model_withoutrandomeffect) 
## [1] 4007.501
#Best model: without random effects 

summary(model_withoutrandomeffect)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *  
##     Genetic_Diversity + Block + (1 | obs)
##    Data: data_G2_all
## 
##      AIC      BIC   logLik deviance df.resid 
##   4007.5   4053.9  -1991.8   3983.5      341 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.42106 -0.07356  0.01527  0.07992  0.27815 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.4413   0.6643  
## Number of obs: 353, groups:  obs, 353
## 
## Fixed effects:
##                                         Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              9.22506    0.51799  17.809  < 2e-16
## Generation_Eggs                         -2.10047    0.28072  -7.483 7.29e-14
## Genetic_DiversityMedium                 -0.65107    0.81458  -0.799    0.424
## Genetic_DiversityLow                     0.10736    0.73631   0.146    0.884
## gen_square                               0.23715    0.03473   6.828 8.59e-12
## BlockBlock4                             -0.12145    0.08305  -1.462    0.144
## BlockBlock5                             -0.48597    0.09723  -4.998 5.79e-07
## Generation_Eggs:Genetic_DiversityMedium  0.32468    0.44435   0.731    0.465
## Generation_Eggs:Genetic_DiversityLow    -0.24262    0.40166  -0.604    0.546
## Genetic_DiversityMedium:gen_square      -0.06201    0.05501  -1.127    0.260
## Genetic_DiversityLow:gen_square          0.01594    0.04970   0.321    0.748
##                                            
## (Intercept)                             ***
## Generation_Eggs                         ***
## Genetic_DiversityMedium                    
## Genetic_DiversityLow                       
## gen_square                              ***
## BlockBlock4                                
## BlockBlock5                             ***
## Generation_Eggs:Genetic_DiversityMedium    
## Generation_Eggs:Genetic_DiversityLow       
## Genetic_DiversityMedium:gen_square         
## Genetic_DiversityLow:gen_square            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnrt_E Gnt_DM Gnt_DL gn_sqr BlckB4 BlckB5 G_E:G_DM G_E:G_DL
## Genrtn_Eggs -0.972                                                            
## Gntc_DvrstM -0.629  0.618                                                     
## Gntc_DvrstL -0.696  0.683  0.442                                              
## gen_square   0.936 -0.989 -0.595 -0.658                                       
## BlockBlock4 -0.097  0.000  0.011  0.007  0.000                                
## BlockBlock5 -0.099  0.012  0.013  0.023 -0.012  0.466                         
## Gnrt_E:G_DM  0.614 -0.632 -0.978 -0.432  0.625  0.001 -0.007                  
## Gnrt_E:G_DL  0.679 -0.699 -0.432 -0.978  0.691  0.004 -0.003  0.441           
## Gntc_DvrM:_ -0.591  0.625  0.942  0.416 -0.631 -0.001  0.007 -0.989   -0.436  
## Gntc_DvrL:_ -0.654  0.691  0.416  0.942 -0.699 -0.004  0.004 -0.437   -0.989  
##             G_DM:_
## Genrtn_Eggs       
## Gntc_DvrstM       
## Gntc_DvrstL       
## gen_square        
## BlockBlock4       
## BlockBlock5       
## Gnrt_E:G_DM       
## Gnrt_E:G_DL       
## Gntc_DvrM:_       
## Gntc_DvrL:_  0.441
## convergence code: 0
## Model failed to converge with max|grad| = 0.0430026 (tol = 0.002, component 1)
#Note: should be the same as: 
model_lm_withrandomeffects <- lme4::lmer(resid ~ (1|Genetic_Diversity:Generation_Eggs), data = data_G2_all)
lmerTest::rand(model_lm_withrandomeffects)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## resid ~ (1 | Genetic_Diversity:Generation_Eggs)
##                                         npar  logLik    AIC         LRT Df
## <none>                                     3 -67.275 140.55               
## (1 | Genetic_Diversity:Generation_Eggs)    2 -67.275 138.55 -1.4211e-13  1
##                                         Pr(>Chisq)
## <none>                                            
## (1 | Genetic_Diversity:Generation_Eggs)          1
########### 
########### OPTION B: Test with Levene's test on residuals
########### 

# Test variances with Levene's test 
tapply(data_G2_all$resid, list(data_G2_all$Generation_Eggs, data_G2_all$Genetic_Diversity),var)
##          High      Medium         Low
## 2 0.017018250 0.009672648 0.002808944
## 3 0.056168990 0.027477499 0.090848301
## 4 0.096321758 0.171608049 0.178076659
## 5 0.042211123 0.164030698 0.183835255
## 6 0.007039084 0.171966465 0.117733007
# Levene’s test on the varianc
# Using leveneTest()
car::leveneTest(resid ~ interaction(Generation_Eggs, Genetic_Diversity), data = data_G2_all)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group  14  2.4148 0.003093 **
##       338                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# PairedData::levene.var.test(data_G2_all$dat[df$Genetic_Diversity=="A"], data_G2_all$dat[df$Genetic_Diversity=="B"]) 
# 
# library("PairedData") 
# 
# leveneTest(data_G2_all$resid[data_G2_all$Genetic_Diversity=="High"&
#                                     data_G2_all$Generation_Eggs==2],
#                 data_G2_all$resid[data_G2_all$Genetic_Diversity=="Low"&
#                                     data_G2_all$Generation_Eggs==2])
# 
# 
# ?levene.var.test()
# 
# 
# ?median()
# ?leveneTest()
# 

########### 
########### OPTION C: Check heteroscedasticity of residuals 
########### 
#Breusch-Pagan test 

#library(lmtest)
#lmtest::bptest(Initial_model)
#p-value is not significant (i.e., >0.05), we do not reject the null hypothesis. Therefore, we assume that the residuals are homoscedastic

#library(skedastic)
#install.packages("skedastic")
#skedastic::white_lm(Initial_model)





########### 
########### OPTION D: Test with Levene's test on data
########### 

# Test variances with Levene's test 
tapply(data_G2_all$resid, list(data_G2_all$Generation_Eggs, data_G2_all$Genetic_Diversity),var)
##          High      Medium         Low
## 2 0.017018250 0.009672648 0.002808944
## 3 0.056168990 0.027477499 0.090848301
## 4 0.096321758 0.171608049 0.178076659
## 5 0.042211123 0.164030698 0.183835255
## 6 0.007039084 0.171966465 0.117733007
# Levene’s test on the varianc
# Using leveneTest()
car::leveneTest(resid ~ interaction(Generation_Eggs, Genetic_Diversity), data = data_G2_all)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group  14  2.4148 0.003093 **
##       338                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.5.2 Raw data

var_raw_data <- ggplot2::ggplot(data[data$Extinction_finalstatus=="no",], 
                aes(x = factor(Generation_Eggs),  
                                 y = Nb_adults, 
                                 colour = Genetic_Diversity)) + 
  geom_boxplot(outlier.color = NA) +
  geom_point(position = position_dodge(0.75), size = 0.7, alpha = 0.5) +
  scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07", "#E7B800")) +
  ylab("Population size") +
  xlab("Generation") +
  geom_text(x=4, y=450, label="**", size = 5, color = "Black") +
  theme_LO
var_raw_data

# 
# 
# cowplot::save_plot(file = here::here("figures","Variance_rawdata.pdf"),
#                    var_raw_data, base_height = 10/cm(1),           
#                    base_width = 15/cm(1), dpi = 610)






########### 
########### OPTION D: Test with Levene's test on data
########### 

# Test variances with Levene's test 
tapply(data[data$Extinction_finalstatus=="no",]$Nb_adults, 
       list(data[data$Extinction_finalstatus=="no",]$Generation_Eggs,
            data[data$Extinction_finalstatus=="no",]$Genetic_Diversity),var, na.rm=TRUE)
##       High   Medium      Low
## 1    0.000    0.000    0.000
## 2 5442.447 8472.330 5546.986
## 3 6555.977 7844.408 6712.525
## 4 6432.400 4085.585 1243.594
## 5 2658.934 2375.036 1858.627
## 6 1815.000 1940.840 1788.358
# Levene’s test on the varianc
# Using leveneTest()
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no",])
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group  17  9.4055 < 2.2e-16 ***
##       406                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==5,])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.4647 0.6303
##       67
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==4,])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  2  4.9423 0.009952 **
##       67                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==6,])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.3492 0.7065
##       68
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==3,])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.5668   0.57
##       68
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==2,])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  2.0367 0.1383
##       68
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&
                              data$Generation_Eggs==4&data$Genetic_Diversity!="Medium",])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  1  10.689 0.001954 **
##       50                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&
                              data$Generation_Eggs==4&data$Genetic_Diversity!="Low",])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.8783 0.1778
##       42
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity), 
                data = data[data$Extinction_finalstatus=="no"&
                              data$Generation_Eggs==4&data$Genetic_Diversity!="High",])
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  2.1686 0.1483
##       42
tapply(data[data$Extinction_finalstatus=="no"&
              data$Generation_Eggs==4,]$Nb_adults, 
       data[data$Extinction_finalstatus=="no"&
              data$Generation_Eggs==4,]$Genetic_Diversity, 
       var, 
       na.rm=TRUE)
##     High   Medium      Low 
## 6432.400 4085.585 1243.594

3.6 Growth rate G6

3.6.1 Distribution analysis

head(data_phenotyping)
##   Population   Week  Block ID_Rep Divsersity Population.1 Box Initial.N   Start
## 1      High1 5-week Block3     52       High            1   1        30  7/8/21
## 2      High1 5-week Block3     53       High            1   2        30  7/8/21
## 3     High13 5-week Block4    129       High           13   1        30 7/15/21
## 4     High13 5-week Block4    130       High           13   2        30 7/15/21
## 5     High13 5-week Block4    131       High           13   3        30 7/15/21
## 6     High14 5-week Block4    132       High           14   1        30 7/15/21
##       End Larvae Pupae Adults    Lambda Genetic_Diversity Nb_ttx  p_adults
## 1 8/12/21      3     9     97 3.2333333              High    109 0.8899083
## 2 8/12/21      1     1      8 0.2666667              High     10 0.8000000
## 3 8/19/21      2     0     84 2.8000000              High     86 0.9767442
## 4 8/19/21      1     2    104 3.4666667              High    107 0.9719626
## 5 8/19/21      2     1     84 2.8000000              High     87 0.9655172
## 6 8/19/21      1     0     83 2.7666667              High     84 0.9880952
##      p_pupae    p_larvae He_remain  He_start    He_exp    He_end
## 1 0.08256881 0.027522936 0.9986008 0.9986008 0.9679591 0.9666047
## 2 0.10000000 0.100000000 0.9986008 0.9986008 0.9679591 0.9666047
## 3 0.00000000 0.023255814 0.9986008 0.9986008 0.9786327 0.9772634
## 4 0.01869159 0.009345794 0.9986008 0.9986008 0.9786327 0.9772634
## 5 0.01149425 0.022988506 0.9986008 0.9986008 0.9786327 0.9772634
## 6 0.00000000 0.011904762 0.9986008 0.9986008 0.9766045 0.9752381
#Without considering extinct populations 

#tapply(data_phenotyping$Lambda,data_phenotyping$Population,length)

hist(data_phenotyping$Lambda, breaks=50)

# #Which distribution
fitdistrplus::descdist(data_phenotyping$Lambda, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  6.2 
## median:  2.1 
## mean:  2.213082 
## estimated sd:  1.058233 
## estimated skewness:  0.5497908 
## estimated kurtosis:  3.726017
#############################################################
################## Determine distribution ###################
#############################################################
#Test normal
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda, "norm")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
##     1-mle-norm 
## "not rejected"
plot(f1)

f1$aic
## [1] 551.898
#Test log-normal
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda+1, "lnorm")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
##    1-mle-lnorm 
## "not rejected"
plot(f1)

f1$aic
## [1] 552.5385
#Test exp
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda, "exp")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
##  1-mle-exp 
## "rejected"
plot(f1)

f1$aic
## [1] 669.5117
#Test log-normal
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda+1, "gamma")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
##    1-mle-gamma 
## "not rejected"
plot(f1)

f1$aic
## [1] 544.6686
## Poisson lognormal model
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
#summary(mlog)

# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
#summary(msqrt)

# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)

# # Log Normal
# mlogN <- glm(log(Lambda) ~ Genetic_Diversity + Block,  
#                      family = gaussian, data=data_phenotyping)



## Compare AIC
AIC(mlog, msqrt,mnorm)
##       df      AIC
## mlog   7 120.8153
## msqrt  7 156.0901
## mnorm  7 528.8983
#Square root a better fits to the data 
#But we compare different values: Growth rate and Growth rate+1

# ## Simulate data 
x.teo.log <- unlist(simulate(mlog))
x.teo.sqrt <- unlist(simulate(msqrt))
x.teo.norm <- unlist(simulate(mnorm))

# ## QQplot to compared Negative binomial and Poisson log normal distributions
qqplot(data_phenotyping$Lambda, x.teo.log,
       main="QQ-plot", xlab="Observed data", ylab="Simulated data",
       las=1, xlim = c(0, 10), ylim = c(0, 10), col='blue') ## QQ-plot
points(sort(data_phenotyping$Lambda), sort(x.teo.sqrt), pch=1, col='red')
points(sort(data_phenotyping$Lambda), sort(x.teo.norm), pch=1, col='green')
abline(0,1)
# ## Add legend
legend(0, 7, legend=c("Log", "Sqrt", "Normal"),
       col=c("blue", "red", "green"), 
       pch=1, bty="n")

# ## Normal distribution provides a good fit to the data


# QQ plot
# qqnorm(resid(mlog))
# qqline(resid(mlog))
# qqnorm(resid(msqrt))
# qqline(resid(msqrt))
# qqnorm(resid(mnorm))
# qqline(resid(mnorm)) #best fit

#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)

hist(residuals(msqrt),col="yellow",freq=F)  #Seems the better fit

hist(residuals(mnorm),col="yellow",freq=F) #Seems the better fit

#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
##     1-mle-norm 
## "not rejected"
## Compare AIC
AIC(mlog, msqrt, mnorm)
##       df      AIC
## mlog   7 120.8153
## msqrt  7 156.0901
## mnorm  7 528.8983
#SLog a better fits to the data 

3.6.2 Aanalysis

#############################################################
##################        Analysis        ###################
#############################################################

mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)

summary(mod0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##    Data: data_phenotyping
## 
## REML criterion at convergence: 514.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3747 -0.6242 -0.0778  0.5409  3.4451 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.2321   0.4817  
##  Residual               0.7480   0.8649  
## Number of obs: 186, groups:  Population, 57
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              2.569709   0.193717  13.265
## Genetic_DiversityMedium -0.873366   0.238587  -3.661
## Genetic_DiversityLow    -0.700898   0.222333  -3.152
## BlockBlock4              0.220484   0.214487   1.028
## BlockBlock5             -0.003695   0.249586  -0.015
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -0.460                     
## Gntc_DvrstL -0.587  0.363              
## BlockBlock4 -0.636  0.079  0.175       
## BlockBlock5 -0.592  0.089  0.232  0.451
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    5 532.53 548.66 -261.26   522.53                         
## mod0    7 520.90 543.48 -253.45   506.90 15.635  2  0.0004027 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE   df lower.CL upper.CL
##  High                2.64 0.135 42.2     2.31     2.98
##  Medium              1.77 0.198 51.5     1.28     2.26
##  Low                 1.94 0.177 58.2     1.51     2.38
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE   df t.ratio p.value
##  High - Medium    0.873 0.239 48.0  3.653  0.0018 
##  High - Low       0.701 0.223 52.5  3.143  0.0076 
##  Medium - Low    -0.172 0.261 53.0 -0.660  0.7873 
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
############ 
###### CONFIDENCE INTERVAL 
############
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity-1 + Block + (1|Population), 
                   data=data_phenotyping)
summary(mod0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity - 1 + Block + (1 | Population)
##    Data: data_phenotyping
## 
## REML criterion at convergence: 514.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3747 -0.6242 -0.0778  0.5409  3.4451 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.2321   0.4817  
##  Residual               0.7480   0.8649  
## Number of obs: 186, groups:  Population, 57
## 
## Fixed effects:
##                          Estimate Std. Error t value
## Genetic_DiversityHigh    2.569709   0.193717  13.265
## Genetic_DiversityMedium  1.696343   0.227931   7.442
## Genetic_DiversityLow     1.868811   0.190682   9.801
## BlockBlock4              0.220484   0.214487   1.028
## BlockBlock5             -0.003695   0.249586  -0.015
## 
## Correlation of Fixed Effects:
##             Gnt_DH Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM  0.369                     
## Gntc_DvrstL  0.331  0.235              
## BlockBlock4 -0.636 -0.457 -0.441       
## BlockBlock5 -0.592 -0.410 -0.331  0.451
confint(mod0)
##                              2.5 %    97.5 %
## .sig01                   0.1878335 0.6406211
## .sigma                   0.7698401 0.9855240
## Genetic_DiversityHigh    2.2034040 2.9372360
## Genetic_DiversityMedium  1.2630197 2.1273737
## Genetic_DiversityLow     1.5040561 2.2293763
## BlockBlock4             -0.1881822 0.6259567
## BlockBlock5             -0.4815724 0.4683096
############ 
###### PREDICT
############
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_phenotyping$Genetic_Diversity),
                                   Block = "Block4")

data_predict_extinct$predict <- predict(mod0, 
                                        type = "response",
                                        re.form = NA,
                                        newdata = data_predict_extinct)



SUM_GrowthRate <- Rmisc::summarySE(data_phenotyping, 
                        measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)

3.7 Remaining He difference

Rmisc::summarySE(data[data$Generation_Eggs==2,],
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd           se           ci
## 1              High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2            Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3               Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320
mfull <- lm(He_remain ~ Genetic_Diversity, data=data[data$Generation_Eggs==2,])
mless <- lm(He_remain ~ 1, data=data[data$Generation_Eggs==2,])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
##   Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
## 1     81 0.00631                                 
## 2     83 3.14572 -2   -3.1394 20162 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE df lower.CL upper.CL
##  High              0.9918 0.001698 81   0.9877   0.9960
##  Medium            0.6743 0.001840 81   0.6698   0.6788
##  Low               0.5404 0.001513 81   0.5367   0.5441
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.318 0.00250 81 126.829 <.0001 
##  High - Low       0.451 0.00227 81 198.470 <.0001 
##  Medium - Low     0.134 0.00238 81  56.200 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
Rmisc::summarySE(data[data$Generation_Eggs==2,],
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd           se           ci
## 1              High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2            Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3               Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320
#Calcul for end:
Rmisc::summarySE(data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",],
          measurevar = c("He_end"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N    He_end         sd          se          ci
## 1              High 27 0.9663114 0.01931641 0.003717444 0.007641317
## 2            Medium 18 0.6378048 0.03493729 0.008234799 0.017373906
## 3               Low 26 0.5090656 0.03443416 0.006753094 0.013908257
mfull <- lm(He_end ~ Genetic_Diversity, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
mless <- lm(He_end ~ 1, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_end ~ Genetic_Diversity
## Model 2: He_end ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     68 0.06009                                  
## 2     70 2.97522 -2   -2.9151 1649.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE df lower.CL upper.CL
##  High               0.966 0.00572 68    0.952    0.980
##  Medium             0.638 0.00701 68    0.621    0.655
##  Low                0.509 0.00583 68    0.495    0.523
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.329 0.00905 68 36.316  <.0001 
##  High - Low       0.457 0.00817 68 55.978  <.0001 
##  Medium - Low     0.129 0.00912 68 14.124  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.8 Growth rate and He

#Best model with genetic diversity
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)

mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod1 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod2 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod5 <- lme4::lmer(Lambda ~ exp(He_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
MuMIn::model.sel(mod0, mod1,mod2,mod3,mod4,mod5)
## Model selection table 
##       (Int) Blc Gnt_Dvr He_end log(He_end) exp(He_end)             family df
## mod1 0.8154   +          1.766                         gaussian(identity)  6
## mod5 0.3445   +                                 0.8311 gaussian(identity)  6
## mod3 2.5470   +                      1.259             gaussian(identity)  6
## mod0 2.5700   +       +                                gaussian(identity)  7
## mod2 2.0770   +                                        gaussian(identity)  5
## mod4 2.0770   +                                        gaussian(identity)  5
##        logLik  AICc delta weight
## mod1 -257.506 527.5  0.00  0.408
## mod5 -258.097 528.7  1.18  0.226
## mod3 -258.156 528.8  1.30  0.213
## mod0 -257.449 529.5  2.05  0.147
## mod2 -263.551 537.4  9.95  0.003
## mod4 -263.551 537.4  9.95  0.003
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Population'
#Makes sense: best model log He
x <- seq(1:100)
y <- seq(1001:1100)
plot(log(x),y)

#############################################################
##################        Analysis        ###################
#############################################################

mod3 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod4    5 532.53 548.66 -261.26   522.53                         
## mod3    6 521.90 541.26 -254.95   509.90 12.627  1  0.0003802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#############################################################
##################        Predict         ###################
#############################################################

#http://rstudio-pubs-static.s3.amazonaws.com/7024_e3a68a9b35454e74abfe15b621c50502.html
mod3 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
SUM_MOD <- summary(mod3)


stats::confint(mod3)
##                   2.5 %    97.5 %
## .sig01       0.23159556 0.6685194
## .sigma       0.76963964 0.9848593
## (Intercept)  0.09849481 1.5283487
## He_end       0.85203146 2.6836107
## BlockBlock4 -0.19603632 0.6401347
## BlockBlock5 -0.48746506 0.4823031
#confint(mod3, method="Wald")
(estimate_He <- SUM_MOD$coefficients["He_end","Estimate"])
## [1] 1.76556
(IC95_lower <- stats::confint(mod3)["He_end","2.5 %"])
## [1] 0.8520315
(IC95_upper <- stats::confint(mod3)["He_end","97.5 %"])
## [1] 2.683611
#PREDICT
data_predict_lambda <- data.frame(He_end =
  seq(min(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end-0.02),
      max(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end+0.02),
      0.01),
                                   Block = "Block4")



data_predict_lambda$lambda_predict <- predict(mod3, type = "response",
                                        re.form = NA,
                                        newdata = data_predict_lambda)






#Confidence interval
bb <- lme4::bootMer(mod3,
              FUN=function(x)predict(x, 
                                     newdata = data_predict_lambda,
                                     re.form=NA),
              nsim=1000)
# Take quantiles of predictions from bootstrap replicates.
# These represent the confidence interval of the mean at any value of Time.
data_predict_lambda$CI_lower <- apply(bb$t, 2, quantile, 0.025)
data_predict_lambda$CI_upper <- apply(bb$t, 2, quantile, 0.975)



# 
# cbind(data_predict_lambda,predict(mod3, newdata = data_predict_lambda, 
#                                   interval = "confidence", level = 0.95))

# cbind(data_predict_lambda, predict(mod3, newdata = data_predict_lambda, re.form = NA,
#                      interval = "confidence", level = 0.95))
# 


### Plot 
# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_phenotyping,
                            measurevar = c("Lambda"),
                            groupvars = c("Genetic_Diversity",
                                          "He_end",
                                          "Population"), 
                            na.rm=TRUE)

  
PLOT_He_expect <- ggplot2::ggplot(SUM_POP_He_Mean) + 
  ggplot2::geom_ribbon(data = data_predict_lambda, 
                       aes(ymin = CI_lower, 
                           ymax = CI_upper, 
                           x = He_end), 
                       fill="grey80", alpha=0.3) + 
  geom_line(data = data_predict_lambda,
            aes(x = He_end, y = lambda_predict),
            col = "grey40", size = 1.5, linetype = 1) + 
  geom_point(data = SUM_POP_He_Mean, 
             aes(x = He_end,  y = Lambda,
                 colour = Genetic_Diversity,
                 size = N), alpha = 0.8) +
  ylab("Growth rate") +
  xlab("Expected heterozygosity") +
  scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"), 
                    values = c("#00AFBB","#FC4E07","#E7B800")) +
  labs(size="Replicates") + 
  theme_LO
PLOT_He_expect

# 
# #
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"),   PLOT_He_expect,
#           base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)

3.9 Adults G6

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   1818.1   1834.3   -904.1   1808.1      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23572 -0.14532  0.02672  0.14825  0.41394 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.1739   0.4170  
##  Population        (Intercept) 0.1085   0.3294  
## Number of obs: 186, groups:  ID_Rep:Population, 186; Population, 57
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2943     0.0812  52.884  < 2e-16 ***
## Genetic_DiversityMedium  -0.5132     0.1438  -3.569 0.000358 ***
## Genetic_DiversityLow     -0.3047     0.1295  -2.352 0.018660 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.563       
## Gntc_DvrstL -0.625  0.356
# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population/ID_Rep),
#                 family = "poisson", data = data_5week)


#dispersion_lme4::glmer(m0) #dispersion ok 
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1    3 1826.3 1836.0 -910.15   1820.3                        
## m0    5 1818.1 1834.3 -904.07   1808.1 12.164  2   0.002284 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   1818.1   1834.3   -904.1   1808.1      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23572 -0.14532  0.02672  0.14825  0.41394 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.1739   0.4170  
##  Population        (Intercept) 0.1085   0.3294  
## Number of obs: 186, groups:  ID_Rep:Population, 186; Population, 57
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2943     0.0812  52.884  < 2e-16 ***
## Genetic_DiversityMedium  -0.5132     0.1438  -3.569 0.000358 ***
## Genetic_DiversityLow     -0.3047     0.1295  -2.352 0.018660 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.563       
## Gntc_DvrstL -0.625  0.356
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean     SE  df asymp.LCL asymp.UCL
##  High                4.29 0.0812 Inf      4.10      4.49
##  Medium              3.78 0.1188 Inf      3.50      4.06
##  Low                 3.99 0.1011 Inf      3.75      4.23
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.513 0.144 Inf  3.569  0.0010 
##  High - Low       0.305 0.130 Inf  2.352  0.0489 
##  Medium - Low    -0.208 0.156 Inf -1.340  0.3728 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m1    4 1823.6 1836.5 -907.81   1815.6                       
## m0    6 1820.1 1839.5 -904.07   1808.1 7.4783  2    0.02377 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1    3 1844.6 1854.3 -919.30   1838.6                        
## m0    5 1838.5 1854.6 -914.26   1828.5 10.086  2   0.006455 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.10 Adults G2

m0 <- glm(Nb_adults ~ Genetic_Diversity + Block,  family = "poisson", data = data[data$Generation_Eggs==2,])
summary(m0)
## 
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson", 
##     data = data[data$Generation_Eggs == 2, ])
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.7057   -2.7621    0.2126    3.3983   11.4454  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              5.73890    0.01497 383.312  < 2e-16 ***
## Genetic_DiversityMedium -0.19408    0.01628 -11.920  < 2e-16 ***
## Genetic_DiversityLow    -0.19278    0.01460 -13.205  < 2e-16 ***
## BlockBlock4              0.10767    0.01579   6.817  9.3e-12 ***
## BlockBlock5              0.17743    0.01600  11.088  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2379.6  on 83  degrees of freedom
## Residual deviance: 2027.9  on 79  degrees of freedom
## AIC: 2667.2
## 
## Number of Fisher Scoring iterations: 4
m1<- glm(Nb_adults ~ Block,  family = "poisson", data = data[data$Generation_Eggs==2,])
anova(m0,m1,test="Chisq") # difference
## Analysis of Deviance Table
## 
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        79     2027.9                          
## 2        81     2241.4 -2  -213.52 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE  df asymp.LCL asymp.UCL
##  High               5.834 0.01051 Inf     5.809     5.859
##  Medium             5.640 0.01243 Inf     5.610     5.670
##  Low                5.641 0.01022 Inf     5.617     5.666
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate     SE  df z.ratio p.value
##  High - Medium   0.1941 0.0163 Inf 11.920  <.0001 
##  High - Low      0.1928 0.0146 Inf 13.205  <.0001 
##  Medium - Low   -0.0013 0.0161 Inf -0.081  0.9964 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.11 Proportion life stage

##### P_adults 
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 723.41 733.18 -357.70   715.41                     
## mod0    6 722.93 737.58 -355.46   710.93 4.4827  2     0.1063
rm(mod0,mod1)


mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod1    4 893.45 906.31 -442.73   885.45                    
## mod0    6 896.17 915.46 -442.08   884.17 1.284  2     0.5262
rm(mod0,mod1)


##### P_adults 
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 692.89 702.67 -342.45   684.89                     
## mod0    6 692.61 707.26 -340.30   680.61 4.2861  2     0.1173
rm(mod0,mod1)

mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod1    4 539.89 552.75 -265.95   531.89                    
## mod0    6 542.26 561.54 -265.13   530.26 1.639  2     0.4407
rm(mod0,mod1)

##### P_larvae 
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference and convergence issue
## Data: data_phenotyping_4week
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod1    4 441.93 451.71 -216.97   433.93                    
## mod0    6 445.08 459.74 -216.54   433.08 0.853  2     0.6528
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference and convergence issue
## Data: data_phenotyping
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 737.53 750.39 -364.77   729.53                     
## mod0    6 740.62 759.91 -364.31   728.62 0.9082  2      0.635

3.12 Multinomial life stage

# data(alligator)
# head(alligator)
# fit <- multinom(food ~ lake+size+sex, data=alligator, weights=count)
# summary(fit$fitted.values)
data_TEMP <- data_phenotyping_all[,c(1:5,11:13,23)]
data_multinom <- tidyr::gather(data_TEMP,"stage","count",6:8)
data_multinom$Divsersity <- as.factor(data_multinom$Divsersity )
head(data_multinom)
##   Population   Week  Block ID_Rep Divsersity    He_end  stage count
## 1      High1 5-week Block3     52       High 0.9666047 Larvae     3
## 2      High1 5-week Block3     53       High 0.9666047 Larvae     1
## 3     High13 5-week Block4    129       High 0.9772634 Larvae     2
## 4     High13 5-week Block4    130       High 0.9772634 Larvae     1
## 5     High13 5-week Block4    131       High 0.9772634 Larvae     2
## 6     High14 5-week Block4    132       High 0.9752381 Larvae     1
# 
# fit <- nnet::multinom(stage ~ Block + Divsersity + Week + 
#                   ID_Rep , data=data_multinom, weights=count)
# summary(fit$fitted.values)
# fit
# 
# fit0 <- nnet::multinom(stage ~ Block + Divsersity + 
#                    ID_Rep , data=data_multinom, weights=count)
# anova(fit,fit0)

# 
# fit1 <- nnet::multinom(stage ~ Block +  Week + ID_Rep, data=data_multinom, weights=count)
# anova(fit1,fit)
# 
# 
# ### New function to test random effect MCLOGIT
# # data_TEMP <- data_phenotyping_all[,c(1:5,17:19,23)]
# # data_multinom <- tidyr::gather(data_TEMP,"stage","prop",6:8)
# # data_multinom$Divsersity <- as.factor(data_multinom$Divsersity )
# # head(data_multinom)
# 
# 
# #Other example
# library("AER")
# library("mlogit")


#### 1- Dataset 1: Train transport
# data("TravelMode", package="AER")
# head(TravelMode)
# #Here: Each individual can choose across 4 modes of transports: air train bus or car
#   #We can determine if this choice is random 
#     #or determined by different variables as: vcost, travel, etc. 
# 
# TM <- mlogit::mlogit.data(data = TravelMode,   #dataset
#                           choice = "choice",   #Variable containing the choice: yes or no
#                           shape = "long",      #type of  dataset
#                           alt.levels = "mode") #Variable contianing the list of choice 
# 
# fit0 <- mlogit::mlogit.data(data=data_multinom, 
#                             choice = "count", 
#                             alt.levels = "stage", 
#                             shape = "long")
# 
# 




#### 2- Dataset 2: Fishing price 
#https://rdrr.io/cran/mlogit/man/mlogit.html

# data("Fishing", package = "mlogit")
# head(Fishing)
# #There are: 
      #two alternative specific variables : price and catch one individual
      #specific variable (income) 
      #four fishing mode : beach, pier, boat, charter

#A dataframe containing :
  #mode: recreation mode choice, one of : beach, pier, boat and charter,
  #price.beach: price for beach mode
  #price.pier: price for pier mode,
  #price.boat: price for private boat mode,
  #price.charter: price for charter boat mode,
  #catch.beach: catch rate for beach mode,
  #catch.pier: catch rate for pier mode,
  #catch.boat: catch rate for private boat mode,
  #catch.charter: catch rate for charter boat mode,
  #income: monthly income,










#### Discussion with Ann Hess 
## Response is stage
## Weight or group with number 
## As Random effect: ID_Box (6 values)
## As Random effect: Population (one to 6 replicates per population)
## As fixed effect: diversity x Time 
data_multinom$stage <- as.factor(data_multinom$stage)
data_multinom$stage <- factor(data_multinom$stage, levels = c("Larvae", "Pupae", "Adults"))

 
# fit <- nnet::multinom(stage ~ Block + Divsersity + Week + 
#                   ID_Rep , data=data_multinom, weights=count)

# mod0 <- ordinal::clmm(stage ~ Block + Divsersity * Week + 
#                         (1|ID_Rep) + (1|Population), 
#                       data=data_multinom, weights=count)
mod1 <- ordinal::clmm(stage ~ Block + Divsersity + Week +
                        (1|ID_Rep) + (1|Population), 
                      data=data_multinom, weights=count)
mod2 <- ordinal::clmm(stage ~ Block + Divsersity  + 
                        (1|ID_Rep) + (1|Population), 
                      data=data_multinom, weights=count)
mod3 <- ordinal::clmm(stage ~ Block + Week  +
                        (1|ID_Rep) + (1|Population), 
                      data=data_multinom, weights=count)
#anova(mod1,mod0) #sign
anova(mod1,mod2) #sign
## Likelihood ratio tests of cumulative link models:
##  
##      formula:                                                            link:
## mod2 stage ~ Block + Divsersity + (1 | ID_Rep) + (1 | Population)        logit
## mod1 stage ~ Block + Divsersity + Week + (1 | ID_Rep) + (1 | Population) logit
##      threshold:
## mod2 flexible  
## mod1 flexible  
## 
##      no.par   AIC  logLik LR.stat df Pr(>Chisq)    
## mod2      8 16875 -8429.7                          
## mod1      9 14842 -7412.2  2035.1  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod1,mod3) # non sign
## Likelihood ratio tests of cumulative link models:
##  
##      formula:                                                            link:
## mod3 stage ~ Block + Week + (1 | ID_Rep) + (1 | Population)              logit
## mod1 stage ~ Block + Divsersity + Week + (1 | ID_Rep) + (1 | Population) logit
##      threshold:
## mod3 flexible  
## mod1 flexible  
## 
##      no.par   AIC  logLik LR.stat df Pr(>Chisq)
## mod3      7 14838 -7412.2                      
## mod1      9 14842 -7412.2   0.127  2     0.9385

4 SUMMARY ALL ANALYSIS

######### Proba of extinction 
mod0 <- glm(y ~ Genetic_Diversity + Block,  
      data = data_proba_extinction, family = binomial)
mod1 <- glm(y ~ Block ,  
      data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1        81     58.903                        
## 2        79     46.833  2   12.069 0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE  df asymp.LCL asymp.UCL
##  High              -20.07 1855.743 Inf  -4451.10  4410.961
##  Medium             -1.77    0.650 Inf     -3.32    -0.216
##  Low                -1.56    0.532 Inf     -2.83    -0.290
## 
## Results are averaged over the levels of: Block 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate       SE  df z.ratio p.value
##  High - Medium  -18.300 1855.743 Inf -0.010  0.9999 
##  High - Low     -18.506 1855.743 Inf -0.010  0.9999 
##  Medium - Low    -0.206    0.751 Inf -0.274  0.9594 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#########Time of extinction
mod1 <- glm(Generation_Eggs ~ Genetic_Diversity + Block,
                    data = data_timeextinction, family = "poisson")
mod0 <- glm(Generation_Eggs ~  Block,
                    data = data_timeextinction, family = "poisson")
anova(mod0, mod1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        10    0.95459                     
## 2         9    0.74888  1  0.20571   0.6502
#########Dynamic over time 
# Test effect
mod_1 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +   s(Generation_Eggs, k=5) +
                    s(Generation_Eggs, by = Genetic_Diversity, k=5) + 
                    s(Population, bs="re"), method="ML", data=data_dyn)
mod_2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block  +  s(Generation_Eggs, k=5) +
                    s(Population, bs="re"), data=data_dyn, method='ML')

itsadug::compareML(mod_1, mod_2)
## mod_1: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Generation_Eggs, by = Genetic_Diversity, k = 5) + s(Population, 
##     bs = "re")
## 
## mod_2: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) + 
##     s(Population, bs = "re")
## 
## Chi-square test of ML scores
## -----
##   Model    Score Edf Difference    Df p.value Sig.
## 1 mod_2 1975.142   8                              
## 2 mod_1 1974.369  14      0.772 6.000   0.957     
## 
## AIC difference: 2.51, model mod_2 has lower AIC.
#########Growth rate 
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    5 532.53 548.66 -261.26   522.53                         
## mod0    7 520.90 543.48 -253.45   506.90 15.635  2  0.0004027 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE   df lower.CL upper.CL
##  High                2.64 0.135 42.2     2.31     2.98
##  Medium              1.77 0.198 51.5     1.28     2.26
##  Low                 1.94 0.177 58.2     1.51     2.38
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE   df t.ratio p.value
##  High - Medium    0.873 0.239 48.0  3.653  0.0018 
##  High - Low       0.701 0.223 52.5  3.143  0.0076 
##  Medium - Low    -0.172 0.261 53.0 -0.660  0.7873 
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
#########Growth rate vs. He
mod3 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod4    5 532.53 548.66 -261.26   522.53                         
## mod3    6 521.90 541.26 -254.95   509.90 12.627  1  0.0003802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##### P_adults 
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 723.41 733.18 -357.70   715.41                     
## mod0    6 722.93 737.58 -355.46   710.93 4.4827  2     0.1063
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod1    4 893.45 906.31 -442.73   885.45                    
## mod0    6 896.17 915.46 -442.08   884.17 1.284  2     0.5262
##### P_adults 
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 692.89 702.67 -342.45   684.89                     
## mod0    6 692.61 707.26 -340.30   680.61 4.2861  2     0.1173
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod1    4 539.89 552.75 -265.95   531.89                    
## mod0    6 542.26 561.54 -265.13   530.26 1.639  2     0.4407
##### P_larvae 
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod1    4 441.93 451.71 -216.97   433.93                    
## mod0    6 445.08 459.74 -216.54   433.08 0.853  2     0.6528
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 737.53 750.39 -364.77   729.53                     
## mod0    6 740.62 759.91 -364.31   728.62 0.9082  2      0.635

5 All analysis with He start or He

######### Proba of extinction 
#Prepare clean dataset
data_temp <- data[data$Generation_Eggs==6,c("Block","Population","He_start","Extinction_finalstatus")]
data_temp$Extinction_finalstatus <- as.factor(data_temp$Extinction_finalstatus)
# #### HERE ADD 
# data_temp$Extinction_finalstatus <- as.factor(data_temp$Extinction_finalstatus)
# 
# data_temp$y <- ifelse(data_temp$Extinction_finalstatus == "yes", 1, 0)
# 
# 
# mod0 <- glm(y ~ He_start + Block,  data = data_temp, family = binomial)
# mod1 <- glm(y ~ Block ,   data = data_temp, family = binomial)
# anova(mod1, mod0, test = "Chisq")
# summary(mod0)
# 



#########Time of extinction
mod1 <- glm(Generation_Eggs ~ He_start + Block,
                    data = data_timeextinction, family = "poisson")
mod0 <- glm(Generation_Eggs ~  Block,
                    data = data_timeextinction, family = "poisson")
anova(mod0, mod1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ He_start + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        10    0.95459                     
## 2         9    0.76852  1  0.18607   0.6662
summary(mod1)
## 
## Call:
## glm(formula = Generation_Eggs ~ He_start + Block, family = "poisson", 
##     data = data_timeextinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.42430  -0.14557   0.03684   0.07219   0.50431  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   2.4209     1.5248   1.588    0.112
## He_start     -0.9243     2.1583  -0.428    0.668
## BlockBlock4  -0.1473     0.5277  -0.279    0.780
## BlockBlock5  -0.3318     0.4810  -0.690    0.490
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.40752  on 12  degrees of freedom
## Residual deviance: 0.76852  on  9  degrees of freedom
## AIC: 53.687
## 
## Number of Fisher Scoring iterations: 4
#########Dynamic over time 
# Test effect
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*He_start + gen_square*He_start + 
                Block  + (1|obs) + (1|Population), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + He_start + 
                Block  + (1|obs) + (1|Population), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")


#Do not converge
#anova(mod1, mod2, test ="Chisq")  


#########Growth rate 
mod0 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    5 532.53 548.66 -261.26   522.53                         
## mod0    6 521.90 541.26 -254.95   509.90 12.627  1  0.0003802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#########Growth rate vs. He
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_end) + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod4    5 532.53 548.66 -261.26   522.53                         
## mod3    6 522.65 542.01 -255.33   510.65 11.878  1   0.000568 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##### P_adults 
mod0 <- lme4::glmer(p_adults ~ He_end + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod1    4 723.41 733.18 -357.70   715.41                       
## mod0    5 721.94 734.16 -355.97   711.94 3.4655  1    0.06266 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lme4::glmer(p_adults ~ He_end + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 893.45 906.31 -442.73   885.45                     
## mod0    5 894.05 910.13 -442.03   884.05 1.4012  1     0.2365
##### P_adults 
mod0 <- lme4::glmer(p_pupae ~ He_end + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod1    4 692.89 702.67 -342.45   684.89                       
## mod0    5 691.48 703.70 -340.74   681.48 3.4111  1    0.06476 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lme4::glmer(p_pupae ~ He_end + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 539.89 552.75 -265.95   531.89                     
## mod0    5 540.28 556.35 -265.14   530.28 1.6188  1     0.2033
##### P_larvae 
mod0 <- lme4::glmer(p_larvae ~ He_end + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping_4week,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 441.93 451.71 -216.97   433.93                     
## mod0    5 443.34 455.55 -216.67   433.34 0.5929  1     0.4413
mod0 <- lme4::glmer(p_larvae ~ He_end + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"),
                   weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population), 
                   data = data_phenotyping,family = binomial(link = "logit"), 
                    weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ He_end + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    4 737.53 750.39 -364.77   729.53                     
## mod0    5 738.90 754.97 -364.45   728.90 0.6313  1     0.4269